Rev. Mod. Phys. 94, 015005 (2022) - Collective nonlinear dynamics and self-organization in decentralized power grids

I. Introduction

The mitigation of climate change is one of the greatest challenges of mankind. Currently, 65% of all greenhouse gas emissions are caused by the carbon dioxide ( $CO 2$ ${CO}_{2}$ ) emissions from fossil fuel combustion and industrial processes (IPCC, 2014). Thus, a fundamental transformation of our energy system is inevitable to meet the goals of the Paris agreement to limit global warming (Rogelj et al., 2015). Rapid actions are needed as greenhouse gas emissions continue to increase and the remaining carbon budgets shrink at an alarming rate (Figueres et al., 2017; Rockström et al., 2017).

Power grids are at the heart of this transformation. Renewable energy sources have shown remarkable development in recent decades (Edenhofer, Pichs-Madruga, and Sokona, 2011; Wiser et al., 2016; Creutzig et al., 2017), but they fundamentally change the operation of the grids they are connected to. The decarbonization of heating, transport, and other sectors introduces new consumers to the electric power grid. Gas and oil heating is replaced by heat pumps, and district heating plants use resistive heaters in times of high renewable power generation (Gröger, Gasteiger, and Suchsland, 2015). Batteries of electric vehicles have to be charged, which may cause congestion in the grid (Carvalho et al., 2015; Bloess, Schill, and Zerrahn, 2018) but also offer a chance of improving stability (Liu et al., 2013; Gajduk et al., 2014).

Power grids are central to our economy and daily life. It is a “uniquely critical infrastructure,” as a variety of other infrastructures and sectors are dependent on the secure supply of electric power (Van der Vleuten and Lagendijk, 2010). Communication, information processing, public transport, and even cooling infrastructure for food supply cannot work without electric power. The energy transition further increases this dependence, as other sectors are being electrified. Ensuring a stable operation of the power system is thus a necessity and a grand challenge.

Transforming the electric power system to meet the 21st century challenges constitutes a transdisciplinary challenge (Brummitt et al., 2013). Modeling, simulation, and analysis of the physical grid have always been an integral part of power engineering. By now, thousands of simulations are run every hour to assess the state of the grid and potential contingencies. But the grid does not run in isolation: it is affected by the weather, by automatic controls and information systems, by energy markets, and finally also by human operators and consumers. Understanding such interdependent systems is challenging, particularly when they leave the normal state of operation in the case of contingencies. Complexity science can effectively complement detailed simulation models. It can elucidate fundamental mechanisms and interactions, provide explicative models, and identify new and unexpected risks. At the same time large amounts of data become available from new measurement devices and smart grid infrastructures. This fosters an empirical and statistical view on power grid operation and stability. In summary, power systems are becoming increasingly complex and methods from statistical physics, nonlinear dynamics, and network science can play an important part in addressing current challenges.

We now name three central challenges for power system operation arising during the energy transition:

(i) The variability of renewable energy sources is perhaps the most well-known one. Wind and solar power generation is determined by the weather and fluctuates on all temporal and spatial scales. Synoptic and seasonal variability are the most obvious features (Heide et al., 2010; Bloomfield et al., 2018; Staffell and Pfenninger, 2018), but power generation can also fluctuate within seconds due to atmospheric turbulence and the dynamics of clouds (Milan, Wächter, and Peinke, 2013; Anvari, Wächter, and Peinke, 2016; Anvari et al., 2016; Zhang et al., 2019). Fingerprints of these fluctuations can be observed in power grid frequency data (Haehne et al., 2018); cf. also Gorjão, Jumar et al. (2020) for a comprehensive data source.

Longer timescales have attracted increasing interest recently, ranging from interannual (Collins et al., 2018) to decadal variability (Wohland et al., 2019) to the impact of climate change (Wohland et al., 2017; Schlott et al., 2018; Weber et al., 2018). A thorough statistical understanding of all modes of variability and their consequences is integral to system operation.

(ii) Furthermore, wind and solar power are typically generated at locations with favorable natural resources. These locations are often far from consumers, so a long-distance transmission of electric power becomes necessary (Pesch, Allelein, and Hake, 2014). This increases grid loads and makes it vulnerable to large-scale blackouts via cascading failures. Robustness and vulnerability are central topics in network science. Physicists have developed and analyzed a variety of models for cascading failures to understand their propagation, their statistical features, and methods to mitigate them; see Albert, Jeong, and Barabási (2000), Carreras et al. (2002), Motter and Lai (2002), Dobson et al. (2007), Witthaut et al. (2016), Yang, Nishikawa, and Motter (2017), Nesti, Zocca, and Zwart (2018), and Schäfer, Witthaut et al. (2018). More recently the optimal design of networks has attracted strong interest in the scientific community; see Nishikawa and Motter (2006), Katifori, Szöllősi, and Magnasco (2010), and Kaiser, Ronellenfitsch, and Witthaut (2020).

(iii) Finally, renewable power sources are different from conventional power plants operating synchronous machines. They are typically connected to the grid via power electronic inverters with different characteristics and dynamics (Carrasco et al., 2006). Inverters have some important disadvantages, in particular, they have no intrinsic inertia to stabilize grid dynamics (Milano et al., 2018) but offer great flexibility in design and control. Hence, the dynamical stability of complex networked systems is a topic of rapidly growing importance; cf. Anvari, Hellmann, and Zhang (2020). Research must address the stability of existing systems as well as the design of future systems.

The goal of this review is twofold. It first provides a starting point for physicists interested in aspects of power system dynamics, operation, and robustness. To this end, the first part of the review is written as a tutorial. We provide a short review of the main tools from dynamical systems and network science in Sec. II. We then review basic principles of power grid operation and provide an overview of static and dynamic models for power grids in Sec. III. We focus on the mathematical description of the elements of power systems, avoiding most technical details covered in the engineering literature (Grainger and Stevenson, 1994; Kundur, 1994; Machowski, Bialek, and Bumby, 2008; Wood, Wollenberg, and Sheblé, 2014). Section IV then discusses the availability of power system datasets required for own studies and simulations. The tutorial part closes with a review of fundamental aspects of power grid stability in Sec. V.

The second part of the review addresses some recent results obtained at the interface of statistical physics, dynamical systems, network science, and power engineering. We first consider aspects of dynamical stability in Sec. VI. When do stable fixed points exist and how do they depend on local dynamics and network topologies? We then analyze the robustness and vulnerability of power grids in Sec. VII. How does a grid react to the failure of single elements and how do large-scale blackouts emerge? We investigate dynamical perturbations or fluctuations and the transient dynamics of grids in Sec. VIII. We conclude with an outlook on current and future research topics in Sec. IX.

II. Fundamentals of Graphs and Networks

A broad variety of collective dynamical phenomena in future-compliant power systems crucially depend on how consumers, producers, and storage infrastructures are interconnected. The topology of the resulting interactions typically exhibits a variety of complex structural features beyond regular lattices or random structures. We introduce here basic concepts from graph theory used to quantitatively describe such topologies and remark by example where they emerge in the analysis and modeling of power systems. A more comprehensive introduction to graph theory and its applications can be found in a variety of textbooks, such as Bollobás (1998).

A graph $G = (V, E)$ $G = (V, E)$ is given by a set $V = {1, \dots, N}$ $V = {1, \dots, N}$ of vertices together with a set $E$ $E$ of edges that is a symmetric subset of $V \times V$ $V \times V$ ; i.e., if $i, j \in V$ $i, j \in V$ then $(i, j) \in E$ $(i, j) \in E$ if and only if $(j, i) \in E$ $(j, i) \in E$ . Because of this symmetry, such a graph is also called an undirected graph. In a directed graph, the edges are ordered 2-tuples, denoted $(i, j) \in E$ $(i, j) \in E$ , thus have a direction from $j$ $j$ to $i$ $i$ and the edge set is not necessarily symmetric.

If $(i, j) \in E$ $(i, j) \in E$ , we call $i$ $i$ and $j$ $j$ adjacent or neighbors and write $i \sim j$ $i \sim j$ , where for directed graphs $i$ $i$ might be adjacent to $j$ $j$ but not vice versa. The degree of a vertex is the number of its neighbors; for directed graphs the in-degree is the number of edges pointing to the vertex, and the out-degree is the number of edges pointing from the vertex. If the degree of every vertex is the same, $k i \equiv k$ $k_{i} \equiv k$ for all $i$ $i$ , the graph is called $k$ $k$ -regular or simply regular.

A walk in a graph is a sequence of vertices and edges

p = (v 0, e 0, v 1, e 1, \dots, e k - 1, v k),

$p = (v_{0}, e_{0}, v_{1}, e_{1}, \dots, e_{k - 1}, v_{k}),$

(1)

where each successive pair of vertices $v j, v j + 1$ $v_{j}, v_{j + 1}$ is connected by an edge $e j = (v j, v j + 1)$ $e_{j} = (v_{j}, v_{j + 1})$ . If no vertex repeats in such a sequence, it is called a path. A cycle is a walk with at least one edge, with first and last vertices identical and with no other vertices repeating. A graph is connected if there is a path from any of its vertices to any other. A directed graph is strongly connected if there is a directed path from any of its vertices to any other; it is weakly connected if for any two of its vertices $v$ $v$ and $w$ $w$ there is a directed path from $v$ $v$ to $w$ $w$ or from $w$ $w$ to $v$ $v$ . A graph with no cycles is called a forest; a connected forest is a tree.

Paths can be used to unambiguously define a metric on a graph. The unweighted length of a path $p$ $p$ is simply the number $k$ $k$ of edges in the sense of Eq. (1). The geodesic distance between two vertices $v n$ $v_{n}$ and $v m$ $v_{m}$ is defined as the length of the shortest path connecting the two vertices. Note that the shortest path is not necessarily unique; if no such path exists, the distance is infinite. In many applications, such a simple definition of distance is not sufficient, as edges may exhibit heterogeneous features. If a distance or weight $d (e)$ $d (e)$ is assigned to every edge $e \in E$ $e \in E$ , the weighted length of a path is given by

length (p) = k - 1 \sum j = 0 d (e j) .

$length (p) = \sum_{j = 0}^{k - 1} d (e_{j}) .$

(2)

As before the length defines a weighted geodesic distance via the shortest path, provided that all $d (e)$ $d (e)$ are non-negative.

The modeling of power grids and other real-world networks includes but goes far beyond the basic concepts of graph theory. In particular, we may assign quantities such as voltages or power injections to the nodes, currents, and power flows to the edges and study dynamical systems on networks. We will thus need an algebraic description of networks, which links properties of nodes and edges to the previously defined basic properties of graphs. For a more detailed overview of algebraic graph theory and its applications to power grids, see Dörfler, Simpson-Porco, and Bullo (2018).

A power grid has $N$ $N$ vertices, often referred to as buses in power engineering, and $L$ $L$ edges or branches. Properties of vertices such as voltages are encoded in vectors in $R N$ $R^{N}$ or $C N$ $C^{N}$ . The absolute value of each complex number encodes the voltage magnitude and the arguments encode its phase (relative to a reference). Similarly, properties of edges such as flows are encoded in vectors in $R L$ $R^{L}$ or $C L$ $C^{L}$ . The topology of such a network or graph is encoded in the adjacency matrix $A \in R N \times N$ $A \in R^{N \times N}$ with the components

A n, m = {1 if vertices n and m are adjacent, 0 otherwise .

$A_{n, m} = {\begin{cases} 1 & if vertices n and m are adjacent, \\ 0 & otherwise . \end{cases}$

(3)

Furthermore, it is useful to define the node-edge incidence matrix $E \in R N \times L$ $E \in R^{N \times L}$ with the components (Newman, 2010)

E n, ℓ = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 1 if line ℓ starts at node n, - 1 if line ℓ ends at node n, 0 otherwise .

$E_{n, ℓ} = {\begin{cases} 1 & if line ℓ starts at node n, \\ - 1 & if line ℓ ends at node n, \\ 0 & otherwise . \end{cases}$

(4)

We note that power grids abstract to undirected graphs: a transmission line can transmit power in either direction. For each edge we fix an orientation to encode the direction of power flow: positive in one direction and negative in the other one. The orientation, the start, and the end of an edge is arbitrary but must be kept fixed.

The node-edge incidence matrix is useful to keep track of network flows. Suppose that the flows over all edges are denoted as a vector $F \in R L$ $F \in R^{L}$ . The total inflow or outflow at all nodes is then given by $E F \in R N$ $E F \in R^{N}$ . The kernel of $E$ $E$ corresponds to cycle flows: For any vector $F \neq 0$ $F \neq 0$ in the kernel, the inflow and outflow vanish at every node such that the flow must be cyclic. We can fix a basis for the kernel by choosing $L - N + 1$ $L - N + 1$ independent fundamental cycles and encoding this basis in the cycle-edge incidence matrix $C \in R L \times (L - N + 1)$ $C \in R^{L \times (L - N + 1)}$ with components

C ℓ, n = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 1 if edge ℓ = (i, j) belongs to cycle c, - 1 if the reverse edge (j, i) belongs to cycle c, 0 otherwise,

$C_{ℓ, n} = {\begin{cases} 1 & if edge ℓ = (i, j) belongs to cycle c, \\ - 1 & if the reverse edge (j, i) belongs to cycle c, \\ 0 & otherwise, \end{cases}$

(5)

such that $E C = 0$ $E C = 0$ . In a plane graph, a graph drawn on a plane without edge crossings, one typically chooses the facets of the graph as the fundamental cycles. After fixing a basis, any cycle flow (any vector in the kernel of $E$ $E$ ) can be written as $F = C f$ $F = C f$ , with coefficients $f 1, \dots, f L - N + 1$ $f_{1}, \dots, f_{L - N + 1}$ .

Not all edges in a grid are equally strong; for instance, the impedance typically scales with the length of a line such that longer transmission lines will have lower conductances and susceptances. We assign weights $b i j$ $b_{i j}$ to all edges $(i, j)$ $(i, j)$ . Either these weights can be included in a diagonal matrix

B d = diag (b 1, b 2, \dots, b L) \in R L \times L

$B_{d} = diag (b_{1}, b_{2}, \dots, b_{L}) \in R^{L \times L}$

(6)

or we can define a weighted adjacency matrix

A n, m = {b n m if vertices n and m are adjacent, 0 otherwise .

$A_{n, m} = {\begin{cases} b_{n m} & if vertices n and m are adjacent, \\ 0 & otherwise . \end{cases}$

(7)

A particularly important matrix to characterize the stability of a network dynamical system is the Laplacian matrix $Λ \in R N \times N$ $Λ \in R^{N \times N}$ with components (Newman, 2010)

Λ n, m = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ \sum N i = 1 b n i for n = m, - b n m if n \neq m and n, m are adjacent, 0 otherwise .

$Λ_{n, m} = {\begin{cases} \sum_{i = 1}^{N} b_{n i} & for n = m, \\ - b_{n m} & if n \neq m and n, m are adjacent, \\ 0 & otherwise . \end{cases}$

(8)

Using the node-edge incidence matrix, the Laplacian is written as $Λ = E B d E ⊤$ $Λ = E B_{d} E^{⊤}$ , where the superscript $⊤$ $⊤$ denotes the transpose of a matrix. If all weights are positive, $Λ$ $Λ$ is symmetric and positive semidefinite; i.e., all eigenvalues are real and non-negative and can be ordered as

0 = λ 1 \leq λ 2 \leq \dots \leq λ N .

$0 = λ_{1} \leq λ_{2} \leq \dots \leq λ_{N} .$

(9)

The second eigenvalue $λ 2$ $λ_{2}$ measures the algebraic connectivity of a graph (Fiedler, 1973), and the associated eigenvector $v 2$ $v_{2}$ is called the Fiedler vector. More generally, the number of zero eigenvalues of $Λ$ $Λ$ equals the number of connected components of the graph (Newman, 2010).

III. Power System Models

In this section we introduce the basic mathematical models for power system operation and stability. We start with the fundamental relations for currents and power flows in ac power grids leading to static power grid models, then introduce dynamical models for synchronous machines and power electronic inverters.

A. Static models

1. Fundamentals of ac power flow

In an ac power grid the voltages $V$ $V$ and currents $I$ $I$ oscillate approximately sinusoidally with time as follows:

V (t) = V peak cos (ω t + ϑ V) = \sqrt 2 R (V_e i ω t), I (t) = I peak cos (ω t + ϑ I) = \sqrt 2 R (I_e i ω t),

$V (t) = V_{peak} \cos (ω t + ϑ_{V}) = \sqrt{2} ℜ (\underline{V} e^{i ω t}), I (t) = I_{peak} \cos (ω t + ϑ_{I}) = \sqrt{2} ℜ (\underline{I} e^{i ω t}),$

(10)

where we have defined the complex-valued amplitudes as

V_= V peak \sqrt 2 e i ϑ V, I_= I peak \sqrt 2 e i ϑ I .

$\underline{V} = \frac{V_{peak}}{\sqrt{2}} e^{i ϑ_{V}}, \underline{I} = \frac{I_{peak}}{\sqrt{2}} e^{i ϑ_{I}} .$

(11)

Voltage and current are typically out of phase due to capacities or inductances, which is reflected by the phase factors $ϑ V$ $ϑ_{V}$ and $ϑ I$ $ϑ_{I}$ . As the phases of voltage and current play essential roles in describing grid operation, we mainly employ a complex notation in this review. In the engineering literature, the complex quantities $V_$ $\underline{V}$ and $I_$ $\underline{I}$ are referred to as phasors and denoted by an underline. We adopt this notation in the following but denote the imaginary unit by the symbol $i$ $i$ , as is common in physics. The electric power is oscillating too, but only the time-averaged value

P = 1 T \int T 0 V (t) I (t) d t = R (V_I_*),

$P = \frac{1}{T} \int_{0}^{T} V (t) I (t) d t = ℜ (\underline{V} {\underline{I}}^{*}),$

(12)

the real power or active power, can do work. The imaginary part $Q = I (V_I_*)$ $Q = ℑ (\underline{V} {\underline{I}}^{*})$ describes the power temporarily stored in capacitances and inductances and is referred to as the reactive power. Furthermore, one defines the apparent power $S = P + i Q = V_I_*$ $S = P + i Q = \underline{V} {\underline{I}}^{*}$ . These relations show that the relative phases of voltages and currents are essential for the power flow in a grid. We note that deviations from perfect sinusoidal signals do occur in practice: for example, higher harmonics in power electronic sources (Liang and Andalib-Bin-Karim, 2018). Nevertheless, this idealization is extremely helpful for the modeling and analysis of power systems.

The basic relation between voltage at the nodes $k, n \in {1, \dots, N}$ $k, n \in {1, \dots, N}$ of a power grid and the currents flowing between the nodes is given by Ohm’s law:

I_k n = 1 z k n (V_k - V_n) = y k n (V_k - V_n),

${\underline{I}}_{k n} = \frac{1}{z_{k n}} ({\underline{V}}_{k} - {\underline{V}}_{n}) = y_{k n} ({\underline{V}}_{k} - {\underline{V}}_{n}),$

(13)

where $z k n = r k n + i x k n$ $z_{k n} = r_{k n} + i x_{k n}$ is the impedance and $y k n = 1 / z k n$ $y_{k n} = 1 / z_{k n}$ is the admittance of the transmission line between the nodes $k$ $k$ and $n$ $n$ . The admittance is divided into its real and imaginary parts $y k n = g k n + i b k n$ $y_{k n} = g_{k n} + i b_{k n}$ , where $g k n$ $g_{k n}$ is the conductance and $b k n$ $b_{k n}$ is the susceptance. For multicircuit transmission lines, we take $y k n$ $y_{k n}$ to be the sum of the admittances of the single circuits and $y k n = 0$ $y_{k n} = 0$ if no transmission line exists. The total current injected into the node $k$ $k$ is then given by $I_k = \sum n I_k n = \sum n y k n (V_k - V_n)$ ${\underline{I}}_{k} = \sum_{n} {\underline{I}}_{k n} = \sum_{n} y_{k n} ({\underline{V}}_{k} - {\underline{V}}_{n})$ . For actual calculations it is convenient to introduce the nodal admittance matrix $Y \in C N \times N$ $Y \in C^{N \times N}$ with the entries

Y k n = G k n + i B k n = {\sum N n = 1 y k n if k = n, - y k n if k \neq n .

$Y_{k n} = G_{k n} + i B_{k n} = {\begin{cases} \sum_{n = 1}^{N} y_{k n} & if & k = n, \\ - y_{k n} & if & k \neq n . \end{cases}$

(14)

Ohm’s law can then be rewritten in a vectorial form, which yields the network equations

I_= Y V_,

$\underline{I} = Y \underline{V},$

(15)

where $I_= (I_1, I_2, \dots, I_N)$ $\underline{I} = ({\underline{I}}_{1}, {\underline{I}}_{2}, \dots, {\underline{I}}_{N})$ represents the complex currents injected into the $N$ $N$ nodes and $V_= (V_1, V_2, \dots, V_N)$ $\underline{V} = ({\underline{V}}_{1}, {\underline{V}}_{2}, \dots, {\underline{V}}_{N})$ denotes the complex voltages. The apparent power injected to a node $k$ $k$ is then given by

S k = V_k I_* k .

$S_{k} = {\underline{V}}_{k} {\underline{I}}_{k}^{*} .$

(16)

Real-world power transmission and distribution grids are mostly constructed as three-phase systems. Hence, there are three conductors with voltages $V A (t) = V peak cos (ω t + ϑ V)$ $V_{A} (t) = V_{peak} \cos (ω t + ϑ_{V})$ , $V B (t) = V peak cos (ω t + ϑ V + 2 π / 3)$ $V_{B} (t) = V_{peak} \cos (ω t + ϑ_{V} + 2 π / 3)$ , and $V C (t) = V peak cos (ω t + ϑ V + 4 π / 3)$ $V_{C} (t) = V_{peak} \cos (ω t + ϑ_{V} + 4 π / 3)$ to ground. The voltage between two conductors oscillates sinusoidally with an amplitude of $\sqrt 3 V peak$ $\sqrt{3} V_{peak}$ . If all three conductors have the same voltage magnitude and a fixed phase difference, it is sufficient to keep one value of $V_$ $\underline{V}$ and $I_$ $\underline{I}$ . Basic relations such as Eq. (15) remain unchanged, but the transmitted power is 3 times as large as in a single-phase system. The power balance at each node hence reads

S k = 3 V_k I_* k = N \sum n = 1 3 | V k | | V n | (G k n - i B k n) (cos θ k n + i sin θ k n),

$S_{k} = 3 {\underline{V}}_{k} {\underline{I}}_{k}^{*} = \sum_{n = 1}^{N} 3 | V_{k} | | V_{n} | (G_{k n} - i B_{k n}) (\cos θ_{k n} + i \sin θ_{k n}),$

(17)

where $θ k n = θ k - θ n$ $θ_{k n} = θ_{k} - θ_{n}$ is the phase difference of complex voltages $V_k$ ${\underline{V}}_{k}$ and $V_n$ ${\underline{V}}_{n}$ . Real and imaginary parts yield the active and the reactive power balance conditions

P k = N \sum n = 1 3 | V k | | V n | (G k n cos θ k n + B k n sin θ k n),

$P_{k} = \sum_{n = 1}^{N} 3 | V_{k} | | V_{n} | (G_{k n} \cos θ_{k n} + B_{k n} \sin θ_{k n}),$

(18)

Q k = N \sum n = 1 3 | V k | | V n | (G k n sin θ k n - B k n cos θ k n) .

$Q_{k} = \sum_{n = 1}^{N} 3 | V_{k} | | V_{n} | (G_{k n} \sin θ_{k n} - B_{k n} \cos θ_{k n}) .$

(19)

Detailed discussions of the network equations and the resulting power flow equations were provided in Chaps. 7–9 of Grainger and Stevenson (1994) and Chap. 6 of Wood, Wollenberg, and Sheblé (2014).

2. Transformers and the per unit system

Real-world transmission systems are more diverse and include different transmission elements. In particular, we consider the following three types of devices. First, transformers are essential to link different voltage levels, from the transmission grid at highest voltages (380 kV in Europe, up to 765 kV in North America) down to the house connecting lines at low voltages (380 V in Europe). Second, special transformers may shift the phase of voltage and current on the two terminal ends (Verboomen et al., 2005). Third, the model for an ordinary transmission line needs to be extended to take into account the charging capacity typically present in real-world lines.

One can include all these devices in the common network model using a system of rescaled variables, the per unit (pu) system, and the unified model for all transmission elements depicted in Fig. 1(b). The equivalent circuit of this transmission element includes an ideal transformer with the tap ratio $t$ $t$ and the phase shift $θ shift$ $θ^{shift}$ in addition to a $π$ $π$ -equivalent line model. The line model includes the series admittance $y$ $y$ and a charging susceptance $b c$ $b^{c}$ , which is attributed equally to the two end points of the line for simplicity.

Transmission elements. (a) Simple series admittance. (b) Unified transmission element consisting of an ideal transformer and a $π$ $π$ -equivalent line model.

To see how the network equations have to be modified in each case, we first consider a single transmission element as shown in Fig. 1(b). The voltages and currents on the two terminals of the ideal transformer are related by a factor of $N 12 = t 12 e i θ shift 12$ $N_{12} = t_{12} e^{i θ_{12}^{shift}}$ . The current flowing through the series admittance is given by $y 12 (V_2 - V_1 / N 12)$ $y_{12} ({\underline{V}}_{2} - {\underline{V}}_{1} / N_{12})$ according to Ohm’s law. Kirchhoff’s current law for the junctions marked by thick black dots in the figure reads

I_2 = y 21 (V_2 - V_1 N 12) + i 2 b c 12 V_2,

${\underline{I}}_{2} = y_{21} ({\underline{V}}_{2} - \frac{{\underline{V}}_{1}}{N_{12}}) + \frac{i}{2} b_{12}^{c} {\underline{V}}_{2},$

(20)

N * 12 I_1 = - y 12 (V_2 - V_1 N 12) + i 2 b c 12 V_1 N 12 .

$N_{12}^{*} {\underline{I}}_{1} = - y_{12} ({\underline{V}}_{2} - \frac{{\underline{V}}_{1}}{N_{12}}) + \frac{i}{2} b_{12}^{c} \frac{{\underline{V}}_{1}}{N_{12}} .$

(21)

It is convenient to introduce scaled units, which are referred to as the “per unit” or “pu” system in power engineering. We fix a reference value $S base$ $S_{base}$ for the power in the entire grid and a reference value $V base, k$ $V_{base, k}$ for the voltage separately for every voltage level $k$ $k$ . The reference values for the currents and admittances are then given by $I base, k = S base / (3 V base, k)$ $I_{base, k} = S_{base} / (3 V_{base, k})$ and $Y base, k = I base, k / V base, k$ $Y_{base, k} = I_{base, k} / V_{base, k}$ . Generally, one selects the nominal voltage of each voltage level as the reference value, while a typical value for the power would be $S base = 100 MVA$ $S_{base} = 100 MVA$ . For the $π$ $π$ -transmission line depicted in Fig. 1, this choice yields $V base, 1 / V base, 2 = t 12$ $V_{base, 1} / V_{base, 2} = t_{12}$ . Dividing Eq. (20) by $I base, 2 = Y base, 2 V base, 2$ $I_{base, 2} = Y_{base, 2} V_{base, 2}$ then yields

I_2 I base, 2 = y 21 Y base, 2 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ V_2 V base, 2 - e - i θ shift 12 V_1 t 12 V base, 2      = V base, 1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ + i 2 b c 12 Y base, 2 V_2 V base, 2 \Rightarrow ˜ I_2 = ˜ y 21 (˜ V_2 - e - i θ shift 12 ˜ V_1) + i 2 ˜ b c 12 ˜ V_2,

$\frac{{\underline{I}}_{2}}{I_{base, 2}} = \frac{y_{21}}{Y_{base, 2}} (\frac{{\underline{V}}_{2}}{V_{base, 2}} - e^{- i θ_{12}^{shift}} \frac{{\underline{V}}_{1}}{\underset{= V_{base, 1}}{\underset{⏟}{t_{12} V_{base, 2}}}}) + \frac{i}{2} \frac{b_{12}^{c}}{Y_{base, 2}} \frac{{\underline{V}}_{2}}{V_{base, 2}} \Rightarrow {\tilde{\underline{I}}}_{2} = {\tilde{y}}_{21} ({\tilde{\underline{V}}}_{2} - e^{- i θ_{12}^{shift}} {\tilde{\underline{V}}}_{1}) + \frac{i}{2} {\tilde{b}}_{12}^{c} {\tilde{\underline{V}}}_{2},$

(22)

where we use a tilde to denote the scaled quantities in the pu system, as with $˜ I_2 = I_2 / I base, 2$ ${\tilde{\underline{I}}}_{2} = {\underline{I}}_{2} / I_{base, 2}$ . Equation (21) can be rescaled analogously if we adopt the convention that $θ shift 12 = - θ shift 21$ $θ_{12}^{shift} = - θ_{21}^{shift}$ . We note that $Y base, 2$ $Y_{base, 2}$ and not $Y base, 1$ $Y_{base, 1}$ has to be used in the normalization of the impedance $˜ y 12 = ˜ y 21 = y 12 / Y base, 2$ ${\tilde{y}}_{12} = {\tilde{y}}_{21} = y_{12} / Y_{base, 2}$ , as the line is on the 2 side of the transformer, not on the 1 side.

In a large power grid with many nodes and transmission lines, we then obtain

˜ I_k = N \sum n = 1 ˜ y k l (˜ V_k - e i θ shift k n ˜ V_n) + i 2 ˜ b c k n ˜ V_k .

${\tilde{\underline{I}}}_{k} = \sum_{n = 1}^{N} {\tilde{y}}_{k l} ({\tilde{\underline{V}}}_{k} - e^{i θ_{k n}^{shift}} {\tilde{\underline{V}}}_{n}) + \frac{i}{2} {\tilde{b}}_{k n}^{c} {\tilde{\underline{V}}}_{k} .$

(23)

We note that an ideal transformer appears as a simple series admittance $˜ y 12$ ${\tilde{y}}_{12}$ in the pu system. This simplification will be used in the modeling of generators, which are typically connected to the grid via a step-up transformer.

Rewriting Eq. (23) in a vectorial form, we recover the network equations in the form of Eq. (15), however, with a modified nodal admittance matrix $˜ Y$ $\tilde{Y}$ defined as

˜ Y k n = {\sum N n = 1 ˜ y k n + i ˜ b c k n / 2 if k = n, - ˜ y k n e i θ shift k l if k \neq n .

${\tilde{Y}}_{k n} = {\begin{cases} \sum_{n = 1}^{N} {\tilde{y}}_{k n} + i {\tilde{b}}_{k n}^{c} / 2 & if & k = n, \\ - {\tilde{y}}_{k n} e^{i θ_{k l}^{shift}} & if & k \neq n . \end{cases}$

(24)

Finally, the power is rescaled to the pu system by dividing Eq. (17) by the reference value $S base = 3 V base, k I base, k$ $S_{base} = 3 V_{base, k} I_{base, k}$ , which yields $˜ S k = ˜ V_k ˜ I_* k$ ${\tilde{S}}_{k} = {\tilde{\underline{V}}}_{k} {\tilde{\underline{I}}}_{k}^{*}$ as desired.

The pu system thus simplifies the network equations, as we can treat all transmission elements alike. In addition, it is advantageous for numerical calculations, as all quantities are of order unity such that rounding errors are minimized. All details of the transmission elements are absorbed into the nodal admittance matrix $˜ Y$ $\tilde{Y}$ . In the absence of phase shifters, $˜ Y = ˜ Y ⊤$ $\tilde{Y} = {\tilde{Y}}^{⊤}$ is symmetric. A detailed discussion of the pu systems and their benefits was provided in Chaps. 1 and 2 of Grainger and Stevenson (1994).

3. The ac load flow equations

Thus far we have established a set of nonlinear algebraic equations describing the steady state of ac power grids, linking nodal power injections $P n$ $P_{n}$ and $Q n$ $Q_{n}$ to nodal state variables $θ n$ $θ_{n}$ and $| V n |$ $| V_{n} |$ . Two of these four quantities are fixed externally for every node, for instance, via power demand. The remaining unknown variables are obtained by numerical solution of the load flow equations; see Chap. 9 of Grainger and Stevenson (1994) and Chap. 6 of Wood, Wollenberg, and Sheblé (2014) for detailed introductions.

In practice, we distinguish three types of nodes. A $P V$ $P V$ bus is typically connected to a generator, which provides a fixed power output $P$ $P$ at fixed voltage magnitudes $| V |$ $| V |$ . A $P Q$ $P Q$ bus represents a node with a given consumption, such that the net injected active power $P$ $P$ and reactive power $Q$ $Q$ are given, while the voltages at these nodes are unknown. A special kind of bus is the slack bus, which captures the supply of the power necessary to have overall power balance in the system to ensure that there is a steady-state solution. It acts as an ideal voltage source (fixed $V_$ $\underline{V}$ ) where the parameters $P$ $P$ and $Q$ $Q$ remain unspecified to balance active and reactive power during the iteration toward the steady-state solution. This is necessary because power losses on the transmission lines are not known a priori before the solution is obtained.

Suppose that we are interested in a transmission grid of $N ℓ$ $N_{ℓ}$ loads, $N g$ $N_{g}$ generators, plus one generator node taken as slack. We then have $2 N ℓ + N g$ $2 N_{ℓ} + N_{g}$ unknown state variables: the phases of all $N g + N ℓ$ $N_{g} + N_{ℓ}$ nonslack nodes and the voltage magnitudes of $N ℓ$ $N_{ℓ}$ load nodes. These unknown variables are fixed by $2 N ℓ + N g$ $2 N_{ℓ} + N_{g}$ nonlinear algebraic equations as follows:

P n (θ 1, \dots, θ N g + N ℓ, | V 1 |, \dots, | V N ℓ |) = P sp n, Q k (θ 1, \dots, θ N g + N ℓ, | V 1 |, \dots, | V N ℓ |) = Q sp k,

$P_{n} (θ_{1}, \dots, θ_{N_{g} + N_{ℓ}}, | V_{1} |, \dots, | V_{N_{ℓ}} |) = P_{n}^{sp}, Q_{k} (θ_{1}, \dots, θ_{N_{g} + N_{ℓ}}, | V_{1} |, \dots, | V_{N_{ℓ}} |) = Q_{k}^{sp},$

(25)

where $n$ $n$ labels all nonslack nodes, while $k$ $k$ labels load nodes only. The superscript “sp” indicates that values on the right side are specified beforehand and the functions $P n (\cdot)$ $P_{n} (\cdot)$ and $Q k (\cdot)$ $Q_{k} (\cdot)$ are given by Eqs. (18) and (19).

A common and effective way to solve this equation system is the Newton-Raphson method [cf. Chap. 9 of Grainger and Stevenson (1994) and Chap. 6 of Wood, Wollenberg, and Sheblé (2014)], which iteratively updates the state vector from an initial guess at the solution of Eqs. (25). The solution for the state vector describes the power grid in steady-state operation. We note that a set of nonlinear equations may exhibit a single well-defined solution but may also have no or several solutions. The lack of a solution indicates an unstable situation that is further discussed in Sec. V.A, while multistability is addressed in Sec. VI.C.3. Recent advances and challenges in numerical methods for ac power flow computations were discussed by Trias (2012) and Mehta, Molzahn, and Turitsyn (2016).

4. Linearized power flow

For small loads and losses in a power grid, load flow calculations simplify considerably [see Chap. 6.18 of Wood, Wollenberg, and Sheblé (2014)] via linearized power flow equations that are based on three approximations. (i) For each transmission element, the resistance and the charging susceptance are small compared to the reactance and are thus neglected. Hence, the admittance is purely imaginary [ $y n k = 1 / (i x n k)$ $y_{n k} = 1 / (i x_{n k})$ ]. (ii) Variations of the voltage magnitude are typically small in transmission grids, such that we can fix them at the reference value of the respective voltage level. In the pu system we thus write $V_n = e i θ n$ ${\underline{V}}_{n} = e^{i θ_{n}}$ , and the power balance at a node $n$ $n$ reads [cf. Eqs. (18) and (20)]

S_n = V_n I_* n = N \sum k = 1 1 i ˜ x n k (V_* n - e - i θ shift n k V_* k) V_n .

${\underline{S}}_{n} = {\underline{V}}_{n} {\underline{I}}_{n}^{*} = \sum_{k = 1}^{N} \frac{1}{i {\tilde{x}}_{n k}} ({\underline{V}}_{n}^{*} - e^{- i θ_{n k}^{shift}} {\underline{V}}_{k}^{*}) {\underline{V}}_{n} .$

(26)

Technically, all nodes must then be considered $P V$ $P V$ buses, such that no equations for the reactive power must be taken into account. Taking the real part of Eq. (26), we obtain the following balance equation for the real power:

P s p n = N \sum k = 1 P n \to k = N \sum k = 1 1 ˜ x n k sin (θ n - θ k - θ shift n, k) .

$P_{n}^{s p} = \sum_{k = 1}^{N} P_{n \to k} = \sum_{k = 1}^{N} \frac{1}{{\tilde{x}}_{n k}} \sin (θ_{n} - θ_{k} - θ_{n, k}^{shift}) .$

(27)

(iii) Finally, small loads imply that the phase differences across a transmission line are small such that the sine function is approximated to first order. The load flow calculations then reduce to the following set of linear equations:

N \sum k = 1 B n k θ k = P sp n - P shift n,

$\sum_{k = 1}^{N} B_{n k} θ_{k} = P_{n}^{sp} - P_{n}^{shift},$

(28)

where $P shift n = \sum k B n k θ shift n k$ $P_{n}^{shift} = \sum_{k} B_{n k} θ_{n k}^{shift}$ accounts for the effects of the phase shifting transformers and the matrix $B \in R N \times N$ $B \in R^{N \times N}$ is the nodal susceptance matrix expressed in per units but without taking into account potential phase shifts:

B n k = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ \sum m ˜ x - 1 n m if n = k, - ˜ x - 1 n k if n is connected to k, 0 otherwise .

$B_{n k} = {\begin{cases} \sum_{m} {\tilde{x}}_{n m}^{- 1} & if n = k, \\ - {\tilde{x}}_{n k}^{- 1} & if n is connected to k, \\ 0 & otherwise . \end{cases}$

(29)

The simplified Eq. (28) is often referred to as the dc approximation, as it is mathematically equivalent to Kirchhoff’s circuit equation for dc electric circuits. Still, it describes the flow of real power in ac power grids. Linear equations can be solved much faster than the nonlinear load flow equations (25), which is advantageous when the flow must be calculated for many different scenarios of generation and load. Furthermore, the linear approximation avoids the problem of multiple or disappearing solutions or multiple optima in optimization problems. Limits of its applicability were discussed in detail in (Purchala et al. (2005), Van Hertem et al. (2006), and Stott, Jardim, and Alsac (2009).

A word of caution is in order, as the symbol $B$ $B$ is used for several related but different quantities. Transmission lines and ordinary transformers have a positive reactance $˜ x > 0$ $\tilde{x} > 0$ . Hence, the off-diagonal elements of the nodal susceptance matrix (29) are negative. In contrast, the imaginary parts of the nodal admittance matrix elements (14) are positive. Both quantities are denoted by the symbol $B$ $B$ , and one must be careful to not confuse them.

5. Matrix formulation

The linearized power flow equations can be condensed in a highly practical compact matrix notation. Let $N$ $N$ denote the number of nodes and $L$ $L$ denote the number of lines in a grid as before. We define the vectors of power injection as $P = (P 1, \dots, P N) ⊤ \in R N$ $P = (P_{1}, \dots, P_{N})^{⊤} \in R^{N}$ , the vector of voltage phase angles as $θ = (θ 1, \dots, θ N) ⊤ \in R N$ $θ = (θ_{1}, \dots, θ_{N})^{⊤} \in R^{N}$ , and the vector of line flows as $F = (F 1, \dots, F L) ⊤ \in R L$ $F = (F_{1}, \dots, F_{L})^{⊤} \in R^{L}$ . As the flow is directed, we must assign an orientation to each line in the grid that is arbitrary but must be kept fixed. We assume that there are no phase-shifting elements and that power injections are balanced; i.e., the $P n$ $P_{n}$ sum to zero.

The directed real power flow $F ℓ$ $F_{ℓ}$ on a line $ℓ$ $ℓ$ from node $m$ $m$ to node $n$ $n$ is given by $F ℓ = x - 1 ℓ (θ m - θ n)$ $F_{ℓ} = x_{ℓ}^{- 1} (θ_{m} - θ_{n})$ . Using the diagonal matrix $B d$ $B_{d}$ defined in Eq. (6) with weights $b ℓ = x - 1 ℓ$ $b_{ℓ} = x_{ℓ}^{- 1}$ and the node-edge incidence matrix (4), the relation of flows and phase angles are given by

F = B d E ⊤ θ .

$F = B_{d} E^{⊤} θ .$

(30)

The real power balance at each node now reads

P = E F,

$P = E F,$

(31)

stating that the power flowing out of each node must equal the power injected at the node. When one combines Eqs. (30) and (31), the equation for the power injections in terms of the voltage angles is obtained as

P = E B d E ⊤ θ = B θ .

$P = E B_{d} E^{⊤} θ = B θ .$

(32)

When Eq. (32) is taken together with Eq. (30), we thus have a linear relation between line flows $F$ $F$ and nodal power injections $P$ $P$ .

The matrix $B = E B d E ⊤$ $B = E B_{d} E^{⊤}$ with the elements given in Eq. (29) is a weighted Laplacian and thus has a single zero eigenvalue if the network is connected; cf. Sec. II. Thus, it is noninvertible and we solve Eq. (32) for $θ$ $θ$ via the Moore-Penrose pseudoinverse $B *$ $B^{*}$ to obtain the line flows directly as a linear function of the nodal power injections

F = B d E ⊤ B *      ≕ PTDF P .

$F = \underset{≕ PTDF}{\underset{⏟}{B_{d} E^{⊤} B^{*}}} P .$

(33)

The matrix prefactor is called the power transfer distribution factor (PTDF) matrix (Wood, Wollenberg, and Sheblé, 2014), as it describes how power injections are distributed throughout the grid.

The zero eigenvalue of the Laplacian $B$ $B$ represents a phase-shift invariance. Uniformly shifting all nodal phases by the same constant $c$ $c$ , $θ n \to θ n + c$ $θ_{n} \to θ_{n} + c$ , does not affect any power flows. Fixing the phase angle $θ k \equiv 0$ $θ_{k} \equiv 0$ at one node $k$ $k$ (the slack node) removes this degree of freedom. We then remove the node from the analysis by removing the $k$ $k$ th row from the vector $θ$ $θ$ and the $k$ $k$ th row and $k$ $k$ th column from the matrix $B$ $B$ . The resulting matrix is called a grounded Laplacian.

6. Generalized linear approximations

Several generalizations of the dc approximation have been developed, and we comment upon three of them.

(i) First, one can improve on the linear approximation of the sine function by rewriting the governing equations of the dc approximation in two parts. We have

P = E B d ψ,

$P = E B_{d} ψ,$

(34)

where $ψ ℓ$ $ψ_{ℓ}$ denotes the sine of the phase difference along the line $ℓ$ $ℓ$ . In the ordinary dc approximation, one simply replaces the sine with its argument; in vectorial form one thus obtains $ψ = E ⊤ θ$ $ψ = E^{⊤} θ$ . If we do not neglect the sine, we instead have

ψ = sin (E ⊤ θ) .

$ψ = \sin (E^{⊤} θ) .$

(35)

The basic idea of the generalized dc approximation (Dörfler and Bullo, 2013b; Simpson-Porco, 2018) is to first obtain $ψ$ $ψ$ approximately and then solve for $θ$ $θ$ . The general solution to the underdetermined equation (34) reads

ψ = E ⊤ B * P + B - 1 d C f,

$ψ = E^{⊤} B^{*} P + B_{d}^{- 1} C f,$

(36)

where $C$ $C$ is the cycle-edge incidence matrix (5) and $f$ $f$ is a vector of cycle flows. To find the actual value of the cycle flows, we would obtain the correct solution of the nonlinear power flow equations, an approach discussed in Sec. VI.C. Here we consider practical approximate solutions and assume the cycle flows $f$ $f$ to be negligible under normal operating conditions. Equation (35) then yields

θ \approx (E E ⊤) - 1 E arcsin (E ⊤ B * P) .

$θ \approx (E E^{⊤})^{- 1} E \arcsin (E^{⊤} B^{*} P) .$

(37)

(ii) A second generalization is to approximately incorporate Ohmic losses in a linear fashion. Two iterative procedures that introduce losses but keep the assumption of fixed voltage magnitudes were described by Stott, Jardim, and Alsac (2009) and Simpson-Porco (2018).

(iii) A third class of linear models includes reactive power flows and voltage variations by appropriate linearization of the nonlinear ac load flow equations. Different approaches were discussed by Coffrin, Van Hentenryck, and Bent (2012), Zhang et al. (2013), and Yang et al. (2019).

7. Hybrid power grids

High voltage dc (HVDC) transmission lines can transmit bulk power over large distances with lower capital cost and lower losses than standard ac lines (Setreus and Bertling, 2008). They are connected to ac via power electronic converter stations, which offers a high degree of flexibility. The transmitted real power can be controlled by the grid operator. Modern converter stations can provide system services such as reactive power compensation. HVDC lines are included in a load flow study rather easily by representing the two converter stations with two additional $P V$ $P V$ buses. One of the converter nodes draws a power $P f < 0$ $P_{f} < 0$ from the grid at node $f$ $f$ , and the other converter then feeds a power

P t = - P f - P loss = - P f - (l 0 - l 1 P f)

$P_{t} = - P_{f} - P_{loss} = - P_{f} - (l_{0} - l_{1} P_{f})$

(38)

back into the grid at node $t$ $t$ . Losses are typically small such that $l 0$ $l_{0}$ is of the order of a few megawatts and $l 1$ $l_{1}$ is of the order of a few percent.

B. Optimal power flow

A fundamental problem in power engineering and energy economics is to determine the cost-optimal dispatch of generators: Which generators should run to satisfy a given demand at minimum costs [see Chap. 13 of Grainger and Stevenson (1994) and Chap. 8 of Wood, Wollenberg, and Sheblé (2014)]? We address this problem here in terms of the dc approximation, taking into account the real power generation and demand but neglecting losses. Suppose that the total demand or load is given by $P load tot$ $P_{tot}^{load}$ and that there are $N g$ $N_{g}$ generators that can be used to satisfy this demand at minimum costs. The actual generation $P gen m$ $P_{m}^{gen}$ of each generator $m = 1, \dots, N g$ $m = 1, \dots, N_{g}$ is bounded by the technical capacity such that we have the constraints

0 \leq P gen m \leq P max m .

$0 \leq P_{m}^{gen} \leq P_{m}^{\max} .$

(39)

The power balance condition in the grid reads $\sum N g m = 1 P gen m = P load tot$ $\sum_{m = 1}^{N_{g}} P_{m}^{gen} = P_{tot}^{load}$ . We assume for simplicity that the costs are proportional to the output $P gen m$ $P_{m}^{gen}$ for each generator. The total variable costs can then be written as

cost tot = N g \sum m = 1 c m P gen m,

${cost}_{tot} = \sum_{m = 1}^{N_{g}} c_{m} P_{m}^{gen},$

(40)

where $c m$ $c_{m}$ denotes the variable costs for generator $m$ $m$ . If there are no further constraints, the solution with minimum costs is simple: Switch on all generators consecutively according to their variable costs $c m$ $c_{m}$ such that

P gen m = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ P max m for c m < c *, P load tot - \sum' m P max m for c m = c *, 0 for c m > c *,

$P_{m}^{gen} = {\begin{cases} P_{m}^{\max} & for & c_{m} < c^{*}, \\ P_{tot}^{load} - \sum_{m}^{'} P_{m}^{\max} & for & c_{m} = c^{*}, \\ 0 & for & c_{m} > c^{*}, \end{cases}$

(41)

where $c *$ $c^{*}$ is the variable cost of the last generator switched on that is identified with the market electricity price. The primed sum in Eq. (41) runs over all generators with $c m < c *$ $c_{m} < c^{*}$ . This is called the merit order principle.

In real-world applications there are many more necessary conditions and constraints. The real power flow in each transmission line must not exceed the line rating in order not to become overloaded and fail:

| F ℓ | \leq F max ℓ for all lines ℓ .

$| F_{ℓ} | \leq F_{ℓ}^{\max} for all lines ℓ .$

(42)

To satisfy this condition we must consider where power is generated and consumed. In particular, the demand or load $P load i$ $P_{i}^{load}$ must be specified separately for every node $i \in {1, \dots, N}$ $i \in {1, \dots, N}$ of the grid. We can then express all line flows $F ℓ$ $F_{ℓ}$ in terms of the dispatch $P gen m$ $P_{m}^{gen}$ as follows. We first introduce the following generator incidence matrix $E gen$ $E^{gen}$ , which indicates where each generator is connected to the grid:

E gen i, m = {1 if generator m is connected to node i, 0 otherwise .

$E_{i, m}^{gen} = {\begin{cases} 1 & if generator m is connected to node i, \\ 0 & otherwise . \end{cases}$

(43)

The resulting net real power injection is thus given by

P n = N g \sum m = 1 E gen n, m P gen m - P load n

$P_{n} = \sum_{m = 1}^{N_{g}} E_{n, m}^{gen} P_{m}^{gen} - P_{n}^{load}$

(44)

for every node $n$ $n$ . Next, the line flows $F ℓ$ $F_{ℓ}$ can be expressed by power injections via Eq. (33), leading to an affine linear relation of flows $F ℓ$ $F_{ℓ}$ and optimization variables $P gen m$ $P_{m}^{gen}$ . Finally, we arrive at the following important optimization problem, which is commonly referred to as dc or linear optimal power flow (OPF):

min P gen m N g \sum m = 1 c m P gen m such that 0 \leq P gen m \leq P max m, | F ℓ | \leq F max ℓ .

$\min_{P_{m}^{gen}} \sum_{m = 1}^{N_{g}} c_{m} P_{m}^{gen} such that 0 \leq P_{m}^{gen} \leq P_{m}^{\max}, | F_{ℓ} | \leq F_{ℓ}^{\max} .$

(45)

Mathematically, this optimization problem is a linear program that admits an efficient solution. Different formulations and their computational efficiency were discussed by Hörsch, Ronellenfitsch et al. (2018). Equation (45) represents the basic optimal power flow problem, and various extensions have been discussed for real-world applications leading to a vast body of literature (Huneault and Galiana, 1991; Momoh, Adapa, and El-Hawary, 1999; Momoh, El-Hawary, and Adapa, 1999; Frank, Steponavice, and Rebennack, 2012a, 2012b). The linear function (40) is a strong simplification of the real cost function. It can be replaced by a convex function without changing the overall complexity of the optimization problem. In contrast, a nonconvex cost function strongly alters the nature of the problem, making it hard to solve in general. Another approach is to consider different modes of operation of a generator. In the simplest case the generator is either on or off. If a generator is on, it typically has a minimum generation such that an additional constraint $P min m \leq P gen m$ $P_{m}^{\min} \leq P_{m}^{gen}$ arises. To take this switching behavior into account, an additional binary optimization variable is introduced, leading to a mixed-integer linear program.

Variable renewable energy sources may be included in two different ways. If renewable generation is taken to be fixed, the real power load at each node $P load i$ $P_{i}^{load}$ is replaced by the residual load, i.e., the difference between the load and renewable generation. However, renewable generation will exceed the load at least temporarily in future energy systems. To allow for a curtailment of the power generation, wind farms and photovoltaic parks are included as generators with zero variable costs and $P max m$ $P_{m}^{\max}$ is set to the assumed or forecasted power generation. For optimization, we may then choose any value $P gen m$ $P_{m}^{gen}$ between $0$ $0$ (complete curtailment) and $P max m$ $P_{m}^{\max}$ , depending on what most benefits the grid. A systematic way of optimizing the dispatch under uncertainties in production and demand is sketched in Sec. VIII.E, where the system constraints are satisfied only with a specified probability.

A further natural extension is to abandon the assumptions of the dc approximation and start with the full nonlinear ac power flow equations. However, the resulting optimization problem is again nonconvex, which makes it much harder to solve. The development of algorithms for this problem is an active field of research; see Erseghe (2014) and Engelmann et al. (2017).

C. Dynamic models

While load flow calculations describe the steady state of a power grid, a large set of models of different complexity is available to analyze its dynamics and the stability of steady states. Which model to choose depends crucially on the phenomena and the timescales to be investigated (Sauer, Pai, and Chow, 2018). The bulk power generation and demand change on timescales of minutes to hours. Optimization models are routinely used to model how the unit commitment and power flows change from hour to hour. In the following, we focus on the next timescale and address aspects of synchronization, transient stability, and, to a lesser extent, voltage dynamics and load-frequency control. These phenomena typically take place on timescales of $10 - 1 - 101 s$ $10^{- 1} - 10^{1} s$ . Faster phenomena, including subtransient effects and electromagnetic propagation effects as well as the modeling of generator control and protection systems, are not covered in this review. Much more detailed and extensive introductions were given in standard textbooks (Kundur, 1994; Machowski, Bialek, and Bumby, 2008; Sauer, Pai, and Chow, 2018). Overviews of different modeling aspects were given by Gajduk, Todorovski, and Kocarev (2014) and Nishikawa and Motter (2015).

1. The swing equation

The classic swing equation describes the dynamics of the mechanical rotation of a synchronous machine; see Chap. 5.1 of Machowski, Bialek, and Bumby (2008) and Nishikawa and Motter (2015). The basic dynamical variable is the mechanical phase angle, which is identical to the voltage phase angle (hence the term “synchronous machine”). The phase angle $δ$ $δ$ is commonly measured relative to a frame of reference rotating at the reference frequency of the grid $ω R$ $ω_{R}$ . The dynamics is then simply given by Newton’s equation as

J ¨ δ = T mech - T el - D mech ω - (D el + κ) Δ ω,

$J \ddot{δ} = T_{mech} - T_{el} - D_{mech} ω - (D_{el} + κ) Δ ω,$

(46)

where $ω = ˙ δ$ $ω = \dot{δ}$ is the angular frequency, $Δ ω = ω - ω R$ $Δ ω = ω - ω_{R}$ is its deviation from the reference, $J$ $J$ is the moment of inertia of the machine, $T mech$ $T_{mech}$ is the mechanical torque driving the rotating machine, and $T el$ $T_{el}$ is the electromagnetic torque due to the power transferred to the grid. The machine experiences mechanical friction with damping coefficient $D mech$ $D_{mech}$ and an effective damping due to damper windings characterized by coefficient $D el$ $D_{el}$ . In addition, one sometimes includes the effects of a frequency controller with gain constant $κ$ $κ$ , which is discussed in Sec. III.D. Aggregating the damping factors $D = D mech + D el + κ$ $D = D_{mech} + D_{el} + κ$ and noting that $˙ δ = ω$ $\dot{δ} = ω$ , we obtain

J ω R ¨ δ + D ω R ˙ δ = ω R ω [(T mech - D mech ω R) ω      = P mech - T el ω    = P el] .

$J ω_{R} \ddot{δ} + D ω_{R} \dot{δ} = \frac{ω_{R}}{ω} [\underset{= P_{mech}}{\underset{⏟}{(T_{mech} - D_{mech} ω_{R}) ω}} - \underset{= P_{el}}{\underset{⏟}{T_{el} ω}}] .$

(47)

We have thus related the dynamics to the net mechanical input power $P mech$ $P_{mech}$ and the electric power $P el$ $P_{el}$ acting on the machine. The swing equation is used to analyze the stability and small-amplitude dynamics around a steady state; hence, we can approximate $ω R / ω \approx 1$ $ω_{R} / ω \approx 1$ .

As for static flows, calculations are often carried out in dimensionless units, i.e., the pu system. To make the swing equation (47) dimensionless, we divide it by the rated power of the machine $P R$ $P_{R}$ and define the quantities

H = (1 / 2) J ω 2 R P R, ˜ D = D ω R P R .

$H = \frac{(1 / 2) J ω_{R}^{2}}{P_{R}}, \tilde{D} = \frac{D ω_{R}}{P_{R}} .$

(48)

The inertia constant $H$ $H$ measures the ratio of stored kinetic energy and output power of the machine when it operates under normal conditions. Typical values are of the order of $4 - 6 s$ $4 - 6 s$ . We then obtain the swing equation

2 H ω R ¨ δ + ˜ D ˙ δ = ˜ P mech - ˜ P el,

$\frac{2 H}{ω_{R}} \ddot{δ} + \tilde{D} \dot{δ} = {\tilde{P}}_{mech} - {\tilde{P}}_{el},$

(49)

where the power is now expressed in the pu system. The motion of the generator is coupled to the rest of the grid via the exchanged real power $P el$ $P_{el}$ , which we discuss later.

2. Principles of synchronous machines

As a next step toward a full machine model, we have to examine how the mechanical rotation induces voltages and currents. A careful introduction was given in Chap. 3.3 of Machowski, Bialek, and Bumby (2008); we provide here only the essential results. Consider first an unloaded round-rotor generator such as that depicted in Fig. 2. The rotor carries a field coil, where a dc current $I f$ $I_{f}$ generates a magnetic field. The stator consists of three coils called A, B, and C in which the rotating field induces an ac voltage; cf. Sec. III.A.1. The magnetic flux in the A coil with mutual inductance $M f$ $M_{f}$ , $Ψ r A (t) = M f I f cos (ω t)$ $Ψ_{r A} (t) = M_{f} I_{f} \cos (ω t)$ , induces the voltage

E r A (t) = - d Ψ r A d t = ω M f I f sin (ω t),

$E_{r A} (t) = - \frac{d Ψ_{r A}}{d t} = ω M_{f} I_{f} \sin (ω t),$

(50)

commonly referred to as the “air-gap” electromagnetic force (EMF). The EMFs in coils B and C are the same up to a phase shift of $\pm 2 π / 3$ $\pm 2 π / 3$ . If the generator is connected to a load or the grid, a current

I A (t) = I M cos (ω t - φ) = \sqrt 2 R (I_e i ω t)

$I_{A} (t) = I_{M} \cos (ω t - φ) = \sqrt{2} ℜ (\underline{I} e^{i ω t})$

(51)

flows through the coil. It is generally phase shifted with respect to the voltage. Here $I_= I m e - i φ / \sqrt 2$ $\underline{I} = I_{m} e^{- i φ} / \sqrt{2}$ is a complex (“phasor”) quantity. The voltage drop at the coil then reads

E r A (t) = ω M f I f sin (ω t) + ω L A    ≕ X A I M cos (ω t - φ) .

$E_{r A} (t) = ω M_{f} I_{f} \sin (ω t) + \underset{≕ X_{A}}{\underset{⏟}{ω L_{A}}} I_{M} \cos (ω t - φ) .$

(52)

An explicit modeling of synchronous machines, including magnetic fluxes and dc currents, is cumbersome and not always necessary for analyzing large grids with tens to thousands of machines. Often an aggregated description in terms of voltage and current phasors is sufficient. Notable exceptions include the modeling of reactive power limits for voltage collapse, the detailed simulation of systemwide transients, or the design of power system stabilizers. In the aggregated description, one defines the two phasor variables $E_f = ω M f I f / \sqrt 2$ ${\underline{E}}_{f} = ω M_{f} I_{f} / \sqrt{2}$ and $E_r$ ${\underline{E}}_{r}$ describing the air-gap EMF. Equation (52) then reads $E_r = E_f - i X A I_$ ${\underline{E}}_{r} = {\underline{E}}_{f} - i X_{A} \underline{I}$ . We can further take into account imperfections, losses, and leakage of magnetic flux by an additional impedance $R l + i X l$ $R_{l} + i X_{l}$ between the air-gap EMF and the terminal voltage of the generator $V_g$ ${\underline{V}}_{g}$ to obtain

V_g = E_f - i X A I_- (R l + i X l) I_.

${\underline{V}}_{g} = {\underline{E}}_{f} - i X_{A} \underline{I} - (R_{l} + i X_{l}) \underline{I} .$

(53)

Hence, the basic relation between the rotor quantities, the air-gap EMF at the stator coils, and the terminal generator voltage can be represented by a simple equivalent circuit, as depicted in Fig. 3. This description is effective and compatible with electric circuit theory.

Schematic of a round-rotor synchronous machine (left image) and a salient pole machine (right image).

Equivalent circuit diagram of (a) a round-rotor synchronous machine and (b) a salient pole machine.

In a salient pole machine the rotor is asymmetric; see Fig. 2. The asymmetry can be taken into account approximately in an effective model that differs only slightly from the case of a round-rotor machine. We decompose all phasors into two contributions, in and out of phase with the magnetic flux of the rotor. The two components are denoted by the subscripts $d$ $d$ and $q$ $q$ referring to the direct and quadrature axes, respectively. To describe the asymmetry, we replace the effective reactance $X A$ $X_{A}$ with $X A d$ $X_{A d}$ and $X A q$ $X_{A q}$ for the $d$ $d$ and $q$ $q$ components. Hence, Eq. (53) is replaced by the two-component equations

E d = V g d + X q I q + R l I d = 0, E q = V g q - X d I d + R l I q = E f,

$E_{d} = V_{g d} + X_{q} I_{q} + R_{l} I_{d} = 0, E_{q} = V_{g q} - X_{d} I_{d} + R_{l} I_{q} = E_{f},$

with $X q = X A q + X l$ $X_{q} = X_{A q} + X_{l}$ and $X d = X A d + X l$ $X_{d} = X_{A d} + X_{l}$ . Figure 3 shows the equivalent diagram of the synchronous machine. Recombining the components yields the phasors of current and terminal voltage as follows:

V_g = V_g d + V_g q = V g d + i V g q, I_= I_d + I_q = I d + i I q .

${\underline{V}}_{g} = {\underline{V}}_{g d} + {\underline{V}}_{g q} = V_{g d} + i V_{g q}, \underline{I} = {\underline{I}}_{d} + {\underline{I}}_{q} = I_{d} + i I_{q} .$

(54)

Why have we introduced the EMFs $E q$ $E_{q}$ and $E p$ $E_{p}$ in addition to $E f$ $E_{f}$ ? The quantities $E d, q$ $E_{d, q}$ characterize the physical source of the induced voltage (the magnetic flux of the field coil), while $E f$ $E_{f}$ quantifies the current in the field coil. The simple relations $E q = E f$ $E_{q} = E_{f}$ and $E d = 0$ $E_{d} = 0$ hold only during the steady operation of the machine. After a disturbance, the magnetic flux changes and so do the effective EMFs $E q, d$ $E_{q, d}$ . Hence, $E q, d$ $E_{q, d}$ become dynamic state variables for which we introduce the equations of motion in the next section. In contrast, $E f$ $E_{f}$ is a fixed system parameter or, if we include excitation control of the synchronous machine, a control variable.

3. Dynamics of synchronous machines

We now proceed to model the transient dynamics of a synchronous machine. A detailed derivation is available in textbooks on electric machines; see Chap. 11.1 of Machowski, Bialek, and Bumby (2008). Here we quote only the result and provide motivation for the form of the resulting equations.

Recall that a machine in the steady state was modeled by an EMF $E_= E_q + E_d$ $\underline{E} = {\underline{E}}_{q} + {\underline{E}}_{d}$ connected to the terminal ends via the reactances $X q$ $X_{q}$ and $X d$ $X_{d}$ , respectively. A similar approach can be used to describe the transient dynamics of a synchronous machine. The EMFs become time dependent, which is commonly indicated by a dash. Similarly, we must assign a different, transient value to the reactances. We stress that these are only effective quantities, and they yield a highly condensed description of the dynamics. A full dynamical model of a synchronous machine must be based on the magnetic fluxes in the different coils, the voltages, and the currents, and the model’s complexity increases quickly.

We first consider the EMF dynamics in the $q$ $q$ axis, which quantifies the flux in the field coil (more precisely the field flux linkage). At steady state the flux is proportional to the voltage exciting the field coil, which leads to $E q = E f$ $E_{q} = E_{f}$ . In the transient regime two additional terms have to be taken into account when one considers the voltage in the field coil. First, the magnitude of flux changes leads to an inductive term proportional to $˙ E' q$ ${\dot{E}}_{q}^{'}$ . Second, there is a term due to the coupling to the stator. A careful accounting of these terms yields the following equation (Machowski, Bialek, and Bumby, 2008):

E' q + T d o ˙ E' q - (X d - X' d) I d = E f .

$E_{q}^{'} + T_{d o} {\dot{E}}_{q}^{'} - (X_{d} - X_{d}^{'}) I_{d} = E_{f} .$

(55)

The following similar relation can be found for the $d$ $d$ -axis EMF:

E' d + T q o ˙ E' d + (X q - X' q) I q = 0 .

$E_{d}^{'} + T_{q o} {\dot{E}}_{d}^{'} + (X_{q} - X_{q}^{'}) I_{q} = 0 .$

(56)

We note that these two equations cover only the transient dynamics after a disturbance. Further models have been introduced to describe the short “subtransient” dynamics. We omit these models, as they typically yield only a small improvement (Stott, 1979; Weckesser, Jóhannsson, and Østergaard, 2013). Note, however, that higher-order models explicitly include damper winding, which we had to include phenomenologically in the swing equation.

4. Synopsis: A hierarchy of dynamical models

Combining the results from Sec. III.C.3 with the swing equation, we obtain a hierarchy of models; see Chap. 11.1 of Machowski, Bialek, and Bumby (2008). In the two-axis or fourth-order model the generator is described by four state variables, the transient EMFs $E' q$ $E_{q}^{'}$ and $E' d$ $E_{d}^{'}$ , the mechanical phase angle $δ$ $δ$ relative to the grid, and its derivative $ω$ $ω$ , which evolve according to

2 H ω R ˙ ω = - ˜ D ω + ˜ P mech - ˜ P el, ˙ δ = ω, T d o ˙ E' q = E f - E' q + (X d - X' d) I d, T q o ˙ E' d = - E' d + (X q - X' q) I q .

$\frac{2 H}{ω_{R}} \dot{ω} = - \tilde{D} ω + {\tilde{P}}_{mech} - {\tilde{P}}_{el}, \dot{δ} = ω, T_{d o} {\dot{E}}_{q}^{'} = E_{f} - E_{q}^{'} + (X_{d} - X_{d}^{'}) I_{d}, T_{q o} {\dot{E}}_{d}^{'} = - E_{d}^{'} + (X_{q} - X_{q}^{'}) I_{q} .$

(57)

The one-axis or third-order model neglects the dynamics of $E' d$ $E_{d}^{'}$ (Schmietendorf et al., 2014). Often it is set to zero ( $E' d = 0$ $E_{d}^{'} = 0$ ), as in the elementary steady-state model discussed in Sec. III.C.2. Hence, the model reduces to three state variables per generator and the last equation in Eq. (57) is omitted.

The second-order or classical model is widely used in the analysis of power system dynamics. It describes the mechanical motion of the generator at constant EMFs:

2 H ω R ˙ ω = - ˜ D ω + ˜ P mech - ˜ P el, ˙ δ = ω .

$\frac{2 H}{ω_{R}} \dot{ω} = - \tilde{D} ω + {\tilde{P}}_{mech} - {\tilde{P}}_{el}, \dot{δ} = ω .$

(58)

To study the resulting dynamics of the state variables, we have to specify the remaining quantities in the equations of motion. The electric quantities $˜ P el$ ${\tilde{P}}_{el}$ , $I q$ $I_{q}$ , and $I d$ $I_{d}$ are time dependent. In fact, they are related to the state variables of all connected generators via algebraic equations describing the grid, as discussed later. The quantities $E f$ $E_{f}$ and $˜ P mech$ ${\tilde{P}}_{mech}$ are set by the control system of the generator, more precisely, the exciter and the governor (Machowski, Bialek, and Bumby, 2008). If the action of the control system is ignored, they become fixed system parameters. The remaining quantities, in particular, the time constants $T$ $T$ and the reactances $X$ $X$ , are regarded as constant system parameters.

We note that detailed models can also include subtransient effects yielding a sixth-order model (Machowski, Bialek, and Bumby, 2008). Furthermore, dynamical models can be reformulated as adaptive networks (Berner, Yanchuk, and Schöll, 2021) or port Hamiltonian systems (Fiaz et al., 2013; Mehrmann et al., 2018), which enables further insights and generalizations.

5. A single generator coupled to an infinite busbar

First basic insights into power system stability can be obtained by focusing on the dynamics of a single generator; see Chaps. 5.3–5.5 of Machowski, Bialek, and Bumby (2008). To this end, we assume that the generator is connected to a large system and that the influence of the generator on that system is negligible. Such a setup is commonly referred to as an “infinite busbar” with a fixed voltage $V s$ $V_{s}$ . We carry out the calculation in the reference frame of the rotor and recall that the angle between the rotor and the stator frame is denoted by $δ$ $δ$ . Hence, the system voltages at the infinite busbar is written as

V s q = V s cos (δ), V s d = - V s sin (δ) .

$V_{s q} = V_{s} \cos (δ), V_{s d} = - V_{s} \sin (δ) .$

(59)

The generator is connected to the infinite busbar through a set of grid elements that are modeled as series impedances. There may first be internal losses or a leakage in the machine that is described by the impedance $R l + i X l$ $R_{l} + i X_{l}$ . There is then typically a step-up transformer that we describe by the impedance $R T + i X T$ $R_{T} + i X_{T}$ , recalling that the pu system takes care of voltage levels automatically. Finally, there are transmission elements with effective impedance $R S + i X S$ $R_{S} + i X_{S}$ to the system. Typically, resistances can be neglected since Ohmic losses as well as shunt impedances are small. The machine EMFs and the system voltages are then related by

V s d = E' d - x' q I q, V s q = E' q + x' d I d, \Rightarrow I q = E' d + V s sin (δ) x' q, I d = - E' q + V s cos (δ) x' d,

$V_{s d} = E_{d}^{'} - x_{q}^{'} I_{q}, V_{s q} = E_{q}^{'} + x_{d}^{'} I_{d}, \Rightarrow I_{q} = \frac{E_{d}^{'} + V_{s} \sin (δ)}{x_{q}^{'}}, I_{d} = \frac{- E_{q}^{'} + V_{s} \cos (δ)}{x_{d}^{'}},$

(60)

where we introduced the shorthand $x' d = X' d + X l + X T + X s$ $x_{d}^{'} = X_{d}^{'} + X_{l} + X_{T} + X_{s}$ and $x' q = X' q + X l + X T + X s$ $x_{q}^{'} = X_{q}^{'} + X_{l} + X_{T} + X_{s}$ . When one neglects Ohmic losses, the air-gap power acting on the machine equals the real power at the system terminal end

˜ P el = V s q I q + V s d I d = E' d I d + E' q I q + (x' d - x' q) I d I q = E' q V s x' d sin (δ) + E' d V s x' q cos (δ) - V 2 s 2 x' q - x' d x' q x' d sin (2 δ) .

${\tilde{P}}_{el} = V_{s q} I_{q} + V_{s d} I_{d} = E_{d}^{'} I_{d} + E_{q}^{'} I_{q} + (x_{d}^{'} - x_{q}^{'}) I_{d} I_{q} = \frac{E_{q}^{'} V_{s}}{x_{d}^{'}} \sin (δ) + \frac{E_{d}^{'} V_{s}}{x_{q}^{'}} \cos (δ) - \frac{V_{s}^{2}}{2} \frac{x_{q}^{'} - x_{d}^{'}}{x_{q}^{'} x_{d}^{'}} \sin (2 δ) .$

(61)

We thus found an explicit expression for the quantities $˜ P el$ ${\tilde{P}}_{el}$ , $I q$ $I_{q}$ , and $I d$ $I_{d}$ in terms of the generator state variables, closing the equations of motion of Sec. III.C.4.

In the widely used classical model of power grid dynamics, one neglects transient saliency and thereby the asymmetry of the rotor geometry in a salient pole machine, assuming that $x' d \approx x' q = x'$ $x_{d}^{'} \approx x_{q}^{'} = x^{'}$ . This assumption simplifies the calculations, as we do not have to separate the $q$ $q$ and $d$ $d$ axes explicitly but can work with phasors. The relation of voltages and currents simplifies to

V_s = E_' - i x' I_\Rightarrow I_= E_- V_s i x',

${\underline{V}}_{s} = {\underline{E}}^{'} - i x^{'} \underline{I} \Rightarrow \underline{I} = \frac{\underline{E} - {\underline{V}}_{s}}{i x^{'}},$

(62)

where the real parts correspond to the $q$ $q$ components and imaginary parts to the $d$ $d$ components. Writing $E_' = E' e - i ϕ$ ${\underline{E}}^{'} = E^{'} e^{- i ϕ}$ and $V_s = V s e - i δ$ ${\underline{V}}_{s} = V_{s} e^{- i δ}$ the real power simplifies to

˜ P el = R (V_s I_*) = E' V s x' sin (δ - ϕ) .

${\tilde{P}}_{el} = ℜ ({\underline{V}}_{s} {\underline{I}}^{*}) = \frac{E^{'} V_{s}}{x^{'}} \sin (δ - ϕ) .$

(63)

Finally, we recall that in the second-order model the dynamics of the rotor EMFs is fully neglected. Hence, $E'$ $E^{'}$ is set to a constant that is chosen to correspond to the steady-state value when the generator provides the power $P + i Q$ $P + i Q$ . We then have to satisfy the relation

P + i Q = V_s I_* = V 2 s - V_s E_i x',

$P + i Q = {\underline{V}}_{s} {\underline{I}}^{*} = \frac{V_{s}^{2} - {\underline{V}}_{s} \underline{E}}{i x^{'}},$

(64)

which can be solved for $E 2$ $E^{2}$ with the result

E 2 = (P x' V s) 2 + (Q x' V s + V s) 2 .

$E^{2} = {(\frac{P x^{'}}{V_{s}})}^{2} + {(\frac{Q x^{'}}{V_{s}} + V_{s})}^{2} .$

(65)

6. Ohmic loads and the Kron reduction

To understand the operation and stability of extended power grids, including collective effects, we extend the model of Sec. III.C.5 to coupled generators. Recall that we have to specify the currents and the real power in terms of the internal EMFs in order to close the equations of motion summarized in Sec. III.C.4. To this end, we have to model the loads and we have to describe the grid explicitly. We formulate Kirchhoff’s equation in a fixed network frame of reference and work out the transformation to the reference frames of the rotating machines explicitly. For simplicity we neglect transient saliency here, setting $x' d = x' q = x'$ $x_{d}^{'} = x_{q}^{'} = x^{'}$ . The derivation follows the presentation by Nishikawa and Motter (2015); a mathematical treatment was provided by Dörfler and Bullo (2013a).

Assume that the generator or active nodes are labeled as $j \in {1, \dots, N g}$ $j \in {1, \dots, N_{g}}$ and the load or passive nodes are labeled as $j \in {N g + 1, \dots, N g + N l}$ $j \in {N_{g} + 1, \dots, N_{g} + N_{l}}$ . For the generator nodes, we have already discussed the relation between internal EMFs $E_'$ ${\underline{E}}^{'}$ and the voltage at the system terminal end $V_$ $\underline{V}$ in Sec. III.C.5. The current injection from internal to terminal end at a node $j$ $j$ is given by

I_in j = E_' j - V_j i x' j .

${\underline{I}}_{j}^{in} = \frac{{\underline{E}}_{j}^{'} - {\underline{V}}_{j}}{i x_{j}^{'}} .$

(66)

Load nodes with a given power demand $S_j = P j + i Q j$ ${\underline{S}}_{j} = P_{j} + i Q_{j}$ are modeled by a fixed admittance to the ground $y ground j = S_* j / V 20$ $y_{j}^{ground} = {\underline{S}}_{j}^{*} / V_{0}^{2}$ . If the voltage magnitude $V j$ $V_{j}$ equals the reference voltage $V 0$ $V_{0}$ , the load nodes then consume the power $V_j I_in * j = V 0 y ground * j V 0 = P j + i Q j$ ${\underline{V}}_{j} {\underline{I}}_{j}^{in *} = V_{0} y_{j}^{ground *} V_{0} = P_{j} + i Q_{j}$ , as desired.

For each node $j \in {1, \dots, N g + N l}$ $j \in {1, \dots, N_{g} + N_{l}}$ we then evaluate Kirchhoff’s current law: The current inflow must equal the current flowing to the other nodes in the grid:

I_in j = N g + N l \sum k = 1 y k j (V_j - V_k) .

${\underline{I}}_{j}^{in} = \sum_{k = 1}^{N_{g} + N_{l}} y_{k j} ({\underline{V}}_{j} - {\underline{V}}_{k}) .$

(67)

We now collect all these linear relations in vectorial form, defining the vectors and matrices

E_= (E_' 1, \dots, E_' N g) ⊤, I_= (I_in 1, \dots, I_in N g) ⊤, V_g = (V_1, \dots, V_N g) ⊤, V_l = (V_N g + 1, \dots, V_N g + N l) ⊤, Y i = diag [(i x' 1) - 1, \dots, (i x' N g) - 1] Y gr g g = diag (y ground N g + 1, \dots, y ground N g + N l) .

$\underline{E} = ({\underline{E}}_{1}^{'}, \dots, {\underline{E}}_{N_{g}}^{'})^{⊤}, \underline{I} = ({\underline{I}}_{1}^{in}, \dots, {\underline{I}}_{N_{g}}^{in})^{⊤}, {\underline{V}}_{g} = ({\underline{V}}_{1}, \dots, {\underline{V}}_{N_{g}})^{⊤}, {\underline{V}}_{l} = ({\underline{V}}_{N_{g} + 1}, \dots, {\underline{V}}_{N_{g} + N_{l}})^{⊤}, Y_{i} = diag [(i x_{1}^{'})^{- 1}, \dots, (i x_{N_{g}}^{'})^{- 1}] Y_{g g}^{gr} = diag (y_{N_{g} + 1}^{ground}, \dots, y_{N_{g} + N_{l}}^{ground}) .$

(68)

Equation (66) and Kirchhoff’s current laws for the generator and load nodes are thus condensed into the form

⎛ ⎜ ⎜ ⎜ ⎝ Y i - Y i 0 - Y i Y g g + Y i Y g l 0 Y l g Y l l + Y gr l l ⎞ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎝ E_V_g V_l ⎞ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎝ I_00 ⎞ ⎟ ⎟ ⎟ ⎠,

$(\begin{matrix} Y_{i} & - Y_{i} & 0 \\ - Y_{i} & Y_{g g} + Y_{i} & Y_{g l} \\ 0 & Y_{l g} & Y_{l l} + Y_{l l}^{gr} \end{matrix}) (\begin{matrix} \underline{E} \\ {\underline{V}}_{g} \\ {\underline{V}}_{l} \end{matrix}) = (\begin{matrix} \underline{I} \\ 0 \\ 0 \end{matrix}),$

(69)

where $Y g g$ $Y_{g g}$ , $Y g l$ $Y_{g l}$ , $Y l g$ $Y_{l g}$ , and $Y l l$ $Y_{l l}$ are the respective partitions of the nodal admittance matrix (14). We can now gradually eliminate the network voltages, using the third row to eliminate $V_l$ ${\underline{V}}_{l}$ and then the second row to eliminate the $V_g$ ${\underline{V}}_{g}$ . This procedure is referred to as Kron reduction, a mathematical accounting of this procedure was discussed by Dörfler and Bullo (2013a). We finally obtain

Y' (1 + Y - 1 i Y') - 1      ≕ Y eff E_= I_,

$\underset{≕ Y^{eff}}{\underset{⏟}{Y^{'} {(1 + Y_{i}^{- 1} Y^{'})}^{- 1}}} \underline{E} = \underline{I},$

(70)

where we defined the shorthand $Y' = Y g g - Y g l (Y l l + Y gr l l) - 1 Y l g$ $Y^{'} = Y_{g g} - Y_{g l} (Y_{l l} + Y_{l l}^{gr})^{- 1} Y_{l g}$ . Equation (70) has the same structure as Eq. (15), but the reduced susceptibility matrix $Y eff$ $Y^{eff}$ represents an effective, not physical network. In summary, Eq. (70) directly yields the currents and the real power acting on the machines in terms of the internal EMFs as follows:

I_in j = N g \sum k = 1 Y eff j k E_' k, P el j = R (V_j I_in * j) = R (E_' j I_in * j) = N g \sum k = 1 R (E_' j Y eff * j k E'_* k) .

${\underline{I}}_{j}^{in} = \sum_{k = 1}^{N_{g}} Y_{j k}^{eff} {\underline{E}}_{k}^{'}, P_{j}^{el} = ℜ ({\underline{V}}_{j} {\underline{I}}_{j}^{in *}) = ℜ ({\underline{E}}_{j}^{'} {\underline{I}}_{j}^{in *}) = \sum_{k = 1}^{N_{g}} ℜ ({\underline{E}}_{j}^{'} Y_{j k}^{eff *} {\underline{E^{'}}}_{k}^{*}) .$

(71)

These quantities are given in the network frame of reference as previously stated. To evaluate them for the dynamical models summarized in Sec. III.C.4, we transfer back to the rotating frame of reference of each machine via

I_j (network frame) = e i δ j (I q j + i I d j) (machine frame), E_' j (network frame) = e i δ j (E' q j + i E' d j) (machine frame) .

${\underline{I}}_{j} (network frame) = e^{i δ_{j}} (I_{q j} + i I_{d j}) (machine frame), {\underline{E}}_{j}^{'} (network frame) = e^{i δ_{j}} (E_{q j}^{'} + i E_{d j}^{'}) (machine frame) .$

Solving for the currents in the machine frame yields

I d j + i I q j = N g \sum k = 1 Y eff j k (E' q k + i E' d k) e i (δ k - δ j),

$I_{d j} + i I_{q j} = \sum_{k = 1}^{N_{g}} Y_{j k}^{eff} (E_{q k}^{'} + i E_{d k}^{'}) e^{i (δ_{k} - δ_{j})},$

(72)

and the real power acting on the $j$ $j$ th machine reads

P el j = R ⎛ ⎝ N g \sum k = 1 Y eff * j k (E' j k + i E' j k) (E' q k - i E' d k) e i (δ j - δ k) ⎞ ⎠ .

$P_{j}^{el} = ℜ (\sum_{k = 1}^{N_{g}} Y_{j k}^{eff *} (E_{j k}^{'} + i E_{j k}^{'}) (E_{q k}^{'} - i E_{d k}^{'}) e^{i (δ_{j} - δ_{k})}) .$

This expression simplifies considerably for the particularly important case of a second-order model. We then have $E' q j = E j = const$ $E_{q j}^{'} = E_{j} = const$ and $E' d j \equiv 0$ $E_{d j}^{'} \equiv 0$ . Writing $Y eff j k = | Y eff j k | e i (γ j k + π / 2)$ $Y_{j k}^{eff} = | Y_{j k}^{eff} | e^{i (γ_{j k} + π / 2)}$ we obtain

P el j = N g \sum k = 1 | Y eff j k | E j E k sin (δ j - δ k - γ j k) .

$P_{j}^{el} = \sum_{k = 1}^{N_{g}} | Y_{j k}^{eff} | E_{j} E_{k} \sin (δ_{j} - δ_{k} - γ_{j k}) .$

(73)

7. The structure-preserving model

In Sec. III.C.6 all loads were modeled as constant admittances to the ground and subsequently eliminated. Bergen and Hill (1981) introduced an alternative model for power grid frequency dynamics, keeping the load nodes and hence the full network structure. Generator nodes are described as in the swing equations, setting the magnitude of the EMF to a constant. The model of load nodes is based on the observation that the real power consumption typically increases with frequency, which is written in the form

P el j (t) = P (0) j + D j ω R ˙ δ j \Rightarrow D j ω R ˙ δ j = P el j - P (0) j .

$P_{j}^{el} (t) = P_{j}^{(0)} + \frac{D_{j}}{ω_{R}} {\dot{δ}}_{j} \Rightarrow \frac{D_{j}}{ω_{R}} {\dot{δ}}_{j} = P_{j}^{el} - P_{j}^{(0)} .$

(74)

We thus obtain an equation of motion for load nodes that is equivalent to the swing equation with a vanishing inertia $H j$ $H_{j}$ . To close the model, one then has to express the real power injections $P el j$ $P_{j}^{el}$ at the generator and load nodes in terms of the voltage phase angles, recalling that the voltage magnitudes are assumed to be constant. This procedure, which includes the elimination of the generator terminal nodes, was described in detail by Nishikawa and Motter (2015). The great advantage of this approach is that the coupling in the resulting equations of motions reflects the true network structure.

8. Dynamics of power electronic inverters for renewable sources

Renewable power sources are generally not equivalent to synchronous machines. All photovoltaic power sources and virtually all wind turbines are connected to the grid via power electronic inverters. The transition to renewable energy supply thus fundamentally changes the dynamics of power grids. In particular, the decrease of inertia provided by conventional synchronous machines is a major challenge for grid stability (Milano et al., 2018).

One generally distinguishes between two modes of operation of power electronic inverters. A grid-following inverter provides a given amount of electric power adjusting to the voltage and frequency provided by the grid. Hence, the operation relies on other resources capable of providing a stable voltage and frequency. In contrast, grid-forming inverters regulate the voltage and frequency to specific set points similar to a synchronous machine. Mixed modes are also possible in principle, regulating either voltage or frequency and following the other.

The development and analysis of new types of grid-forming inverters and inverter-based grids is an active field of research. We thus introduce one important class, the droop-controlled inverter, following Schiffer et al. (2014). The basic state variables of such an inverter are the EMF magnitude $E j$ $E_{j}$ and the phase angle $δ j$ $δ_{j}$ , where $j$ $j$ labels the different generators in the grid. The control system adjusts these state values to maintain the predefined set values of power and frequency. To this end, the control system measures the real and reactive power exchanged with the grid and compares it to predefined set values. In a simple proportional, or droop control, scheme, the frequency control is proportional to the active power deviation and the voltage control is proportional to the reactive power deviation such that we obtain the equations of motion

˙ δ j = ω d - κ act j (P mes j - P d j),

${\dot{δ}}_{j} = ω^{d} - κ_{j}^{act} (P_{j}^{mes} - P_{j}^{d}),$

(75)

T V j ˙ E j = - E j + E d j - κ rec j (Q mes j - Q d j),

$T_{j}^{V} \dot{E_{j}} = - E_{j} + E_{j}^{d} - κ_{j}^{rec} (Q_{j}^{mes} - Q_{j}^{d}),$

(76)

where the superscript “mes” indicates the measured values of real and reactive power and the superscript $d$ $d$ stands for desired values. Naturally, the desired frequency $ω d$ $ω^{d}$ is unique across the grid, whereas the desired voltages $E d j$ $E_{j}^{d}$ may differ. The parameters $κ act j$ $κ_{j}^{act}$ and $κ rec j$ $κ_{j}^{rec}$ are the droop gains for the frequency and voltage, respectively. The measurements are not instantaneous and are characterized by a low-pass filter such that $T j ˙ P mes j = - P mes j + P el j$ $T_{j} {\dot{P}}_{j}^{mes} = - P_{j}^{mes} + P_{j}^{el}$ and $T j ˙ Q mes j = - Q mes j + Q el j$ $T_{j} {\dot{Q}}_{j}^{mes} = - Q_{j}^{mes} + Q_{j}^{el}$ . The recovery time of voltage dynamics is typically much lower than the time constant $T j$ $T_{j}$ of the low-pass filter. We thus simplify the model by setting $T V j = 0$ $T_{j}^{V} = 0$ and obtain (Schiffer et al., 2014)

˙ δ j = ω j,

${\dot{δ}}_{j} = ω_{j},$

(77)

T j ˙ ω j = - ω j + ω d - κ act j (P el j - P d j),

$T_{j} {\dot{ω}}_{j} = - ω_{j} + ω^{d} - κ_{j}^{act} (P_{j}^{el} - P_{j}^{d}),$

(78)

T j ˙ E j = - E j + E d j - κ rec j (Q el j - Q d j) .

$T_{j} {\dot{E}}_{j} = - E_{j} + E_{j}^{d} - κ_{j}^{rec} (Q_{j}^{el} - Q_{j}^{d}) .$

(79)

The real and reactive power exchanged with the grid ( $P e l j$ $P_{j}^{e l}$ and $Q e l j$ $Q_{j}^{e l}$ ) depend on the state of all elements in the grid. To close the equations of motion, one thus has to specify the network equations. Modeling all loads by constant impedances to the ground as in Sec. III.C.6, we find that

P el j + i Q el j = N g \sum k = 1 E_' j Y eff * j k E_' * k,

$P_{j}^{el} + i Q_{j}^{el} = \sum_{k = 1}^{N_{g}} {\underline{E}}_{j}^{'} Y_{j k}^{eff *} {\underline{E}}_{k}^{' *},$

(80)

where $E_j = E j e i δ j$ ${\underline{E}}_{j} = E_{j} e^{i δ_{j}}$ ; cf. Eq. (71).

9. Aggregated dynamical model

The classical model describes the dynamics of power systems as coupled second-order rotators, with each node representing one synchronous generator. Notably, similar equations of motion emerge on much coarser spatial scales (Zhang et al., 1997; You, Vittal, and Wang, 2004; Filatrella, Nielsen, and Pedersen, 2008; Chow, 2013), which emphasizes the generality of the coupled rotator model.

The dynamical model introduced by Filatrella, Nielsen, and Pedersen (2008) considered the aggregated dynamics of regions labeled by $n \in {1, \dots, N}$ $n \in {1, \dots, N}$ . If the internal coupling is sufficiently strong, we take the phase angle to be constant throughout the region. We again measure the phase and frequency relative to a frame rotating at the reference grid frequency $ω R$ $ω_{R}$ while writing $ϕ n (t) = ω R t + δ n (t)$ $ϕ_{n} (t) = ω_{R} t + δ_{n} (t)$ . The equations of motion for $δ n (t)$ $δ_{n} (t)$ are obtained from energy conservation. Each region stores kinetic energy in the rotation of all synchronous machines as follows:

E kin n = J n 2 (˙ ϕ) 2,

$E_{n}^{kin} = \frac{J_{n}}{2} (\dot{ϕ})^{2},$

(81)

where $J n$ $J_{n}$ is the aggregated moment of inertia. At each node a power $P in n$ $P_{n}^{in}$ is injected by generators or withdrawn by consumers and some energy will be dissipated at a rate $P diss n = η n (˙ ϕ n) 2$ $P_{n}^{diss} = η_{n} ({\dot{ϕ}}_{n})^{2}$ . Assuming a constant voltage magnitude, the real power transmitted from a region $n$ $n$ to a region $m$ $m$ is determined by the phase difference

P trans n \to m = K n m sin (δ n - δ m),

$P_{n \to m}^{trans} = K_{n m} \sin (δ_{n} - δ_{m}),$

(82)

with Ohmic losses neglected. $K n m$ $K_{n m}$ is an aggregated quantity summing over all lines connecting the two regions $n$ $n$ and $m$ $m$ . The energy conservation for region $n$ $n$ then gives

P in n = d E kin n d t + P diss n + \sum m P n \to m .

$P_{n}^{in} = \frac{d E_{n}^{kin}}{d t} + P_{n}^{diss} + \sum_{m} P_{n \to m} .$

(83)

Noting that $| ˙ δ | ≪ ω R$ $| \dot{δ} | ≪ ω_{R}$ in the vicinity of the normal operation, we simplify the expression for the change of the kinetic energy $d E kin n / d t \approx ω R J n ¨ δ n$ $d E_{n}^{kin} / d t \approx ω_{R} J_{n} {\ddot{δ}}_{n}$ and the dissipated power $P diss n \approx η ω 2 R + 2 η n ω R ˙ δ n$ $P_{n}^{diss} \approx η ω_{R}^{2} + 2 η_{n} ω_{R} {\dot{δ}}_{n}$ . Substituting these results into the energy conservation law (83) then directly yields the equations of motion (Rohden, 2012)

ω R J n ¨ δ n + 2 ω R η n ˙ δ n = P eff n - \sum m K n m sin (δ n - δ m),

$ω_{R} J_{n} {\ddot{δ}}_{n} + 2 ω_{R} η_{n} {\dot{δ}}_{n} = P_{n}^{eff} - \sum_{m} K_{n m} \sin (δ_{n} - δ_{m}),$

(84)

where $P eff n = P in n - ω 2 R η n$ $P_{n}^{eff} = P_{n}^{in} - ω_{R}^{2} η_{n}$ . Equation (84) has the same structure as the classical second-order model but describes the dynamics on a coarser level. Equation (84) is based on elementary energetic arguments, but similar results can also be obtained via model reduction (Zhang et al., 1997; You, Vittal, and Wang, 2004; Chow, 2013). This emphasizes the generic nature of coupled rotator models in network science.

D. Load-frequency control

A power grid by itself does not store energy, so the generated power must match the demand and the losses at all times. A hierarchy of control systems exists to maintain this balance, which we introduce starting with the swing equation (47). As before, we use a rotating frame of reference such that $ω i = ˙ δ i$ $ω_{i} = {\dot{δ}}_{i}$ denotes the deviation from the nominal grid frequency $ω R$ $ω_{R}$ ,

J i ω R ˙ ω i + D i ω R ω i = P mech i - P el i .

$J_{i} ω_{R} {\dot{ω}}_{i} + D_{i} ω_{R} ω_{i} = P_{i}^{mech} - P_{i}^{el} .$

(85)

In normal operation, all machines $i \in {1, \dots, N}$ $i \in {1, \dots, N}$ in a grid run in synchrony at the same frequency, a fact that we discuss further in Sec. VI. We thus focus on the bulk frequency and define

¯ ¯ ω (t) = ¯ J - 1 N \sum i = 1 J i ω R ω i (t) with ¯ J = N \sum i = 1 J i ω R .

$\bar{ω} (t) = {\bar{J}}^{- 1} \sum_{i = 1}^{N} J_{i} ω_{R} ω_{i} (t) with \bar{J} = \sum_{i = 1}^{N} J_{i} ω_{R} .$

(86)

In many studies, the damping is assumed to be proportional to the inertia of a machines; cf. the Supplemental Material of Motter et al. (2013). Setting $D i = η J i$ $D_{i} = η J_{i}$ we obtain

¯ J ˙ ¯ ω + η ¯ J ¯ ω = Δ P,

$\bar{J} \dot{\bar{ω}} + η \bar{J} \bar{ω} = Δ P,$

(87)

where $Δ P = \sum i P mech i - P el i$ $Δ P = \sum_{i} P_{i}^{mech} - P_{i}^{el}$ is the power balance in the entire grid. We thus see that the power balance directly drives the dynamics of the bulk: A scarcity of generation leads to a decrease, while an overgeneration leads to an increase of the frequency. As the frequency can be measured easily, it is used to control the generation.

Different control mechanisms are distinguished according to the different timescales that they act upon and the different purposes that they serve; see Fig. 4.

(i) The momentary reserve is provided by the inertia of synchronous machines. The higher the inertia $¯ J$ $\bar{J}$ , the slower the grid frequency reacts to a power imbalance. This reaction rate is referred to as the rate of change of frequency (RoCoF) in the engineering literature. Notably renewable power sources with power electronic inverters have no intrinsic inertia, such that frequency stability is an important challenge in the energy transition (Milano et al., 2018).

(ii) Next, the power imbalance is reduced by increasing or decreasing the generation of specific control power plants or, more recently, battery electric storage systems (Fleer and Stenzel, 2016). According to the guidelines of the European Network of Transmission System Operators for Electricity (ENTSO-E) (UCTE Operations Handbook, 2009), primary control is activated as soon as the frequency leaves a small dead band around $ω R$ $ω_{R}$ and must be fully available within 30 s. Up to the dead band, the power is adapted proportionally to the frequency deviation as follows:

P PRL, i (t) = - κ i ω i,

$P_{PRL, i} (t) = - κ_{i} ω_{i},$

(88)

where the parameter $κ i$ $κ_{i}$ describes the sensitivity of primary control for the grid area $i$ $i$ . Primary control stabilizes the frequency but does not restore the nominal grid frequency, which can be seen by adding primary control to the equation of motion (87). For a constant power imbalance $Δ P$ $Δ P$ and a single area we find the fixed point

¯ ω = Δ P κ + η ¯ J .

$\bar{ω} = \frac{Δ P}{κ + η \bar{J}} .$

(89)

The damping terms $\propto D i$ $\propto D_{i}$ in Eq. (85) affect frequency dynamics similar to the primary control $\propto κ i$ $\propto κ_{i}$ , which is often directly included in the swing equation; cf. Sec. III.C.1.

(iii) Secondary control restores the nominal grid frequency $ω R$ $ω_{R}$ and reduces unscheduled power flows between different grid areas on a timescale of minutes. Control power plants adapt their generation according to the following proportional-integral (PI) law (Böttcher et al., 2020):

P SRL, i (t) = - (κ P G i (t - τ) + κ I \int t - τ - \infty G i (t') d t'),

$P_{SRL, i} (t) = - (κ_{P} G_{i} (t - τ) + κ_{I} \int_{- \infty}^{t - τ} G_{i} (t^{'}) d t^{'}),$

(90)

with $κ P$ $κ_{P}$ and $κ I$ $κ_{I}$ tunable gain factors. The local area-control error $G i$ $G_{i}$ is a measure of the power that is missing in area $i$ $i$ . It is determined by the difference between the expected primary control power and the deviations $Δ F i$ $Δ F_{i}$ of the power flow to neighboring control areas, which is given as $G i = κ i ω i - Δ F i$ $G_{i} = κ_{i} ω_{i} - Δ F_{i}$ . Note that while the PI controller is linear, the local area-control error $G i$ $G_{i}$ depends nonlinearly on the system state. The finite time delay $τ$ $τ$ in the control [Eq. (90)] results from the time required for the measurement of $G i$ $G_{i}$ , communication, and the adaptation of the control power generation. The secondary control requires measuring local area-control errors $G i$ $G_{i}$ with a cycle time of 2–5 s (UCTE Operations Handbook, 2009). The delay time $τ$ $τ$ can therefore be expected to be of similar magnitude but may vary in different control areas and with time (Böttcher et al., 2020).

(iv) Tertiary control power is invoked on even longer timescales. In Europe, the tertiary control power must be activated within at most 15 min. Tertiary control is partly manually operated and used mostly to restore the automatic control reserve.

Load-frequency control scheme according to ENTSO-E (UCTE Operations Handbook, 2009).

IV. Power Grid Topologies and Datasets

We previously introduced the various components of power grids, in particular, power lines and assorted producers and consumers at the nodes. To derive a model of a complete power grid, the missing ingredient is the topology, that is, how these components are connected. Besides studying concrete examples of real-world power grids, it is also desirable to obtain plausible synthetic networks to allow a more systematic study of the impact of topological features on dynamical properties.

A. General aspects

Large power grids, especially at the continental scale, operate at many different voltage levels. Long-range transmission is achieved through high voltage connections in order to minimize losses. In Central Europe, for example, the transmission grid operates at 220 and 380 kV. The regional distribution is then achieved by grids at successively lower voltages that are connected to the high voltage transmission grid at substations (transformers). Typical voltage levels in Europe range from high voltages (such as 60 and 110 kV) over medium voltages (3–30 kV) to low voltages (such as 230 and 380 V). Typical line parameters for transmission grids are summarized in Table I.

TABLE I.

Typical parameters of overhead transmission lines at different voltage levels for the United Kingdom and the USA according to Kundur (1994) and Machowski, Bialek, and Bumby (2008) and for Germany according to Oeding and Oswald (2016). Ohmic resistance, susceptance, and charging capacity are proportional to the length of the line and are given per km.

$f n$ $f_{n}$ (Hz)	$V n$ $V_{n}$ (kV)	$r$ $r$ ( $Ω / km$ $Ω / km$ )	$x$ $x$ ( $Ω / km$ $Ω / km$ )	$b c$ $b_{c}$ ( $μ S / km$ $μ S / km$ )
50 (UK)	400	0.018	0.265	5.36
	230	0.05	0.488	3.371
50 (Germany)	380	0.03	0.246	4.335
	220	0.06	0.301	3.927
60 (USA)	345	0.037	0.367	4.518
	500	0.028	0.325	5.2
	765	0.012	0.329	4.978

Grids at these different voltage levels are constructed in different ways and with different trade-offs and goals in mind and thus exhibit different network topologies. A low voltage distribution grid is typically organized as a tree graph, which in this context is called a radial design. While loops might physically exist, there often is a switch interrupting them. We remark that only if a line fails is the loop closed to resume supply downstream of the failure.

High voltage transmission grids require uninterrupted operation even in cases where a line fails (the $N - 1$ $N - 1$ criterion). This implies that the network is meshed; i.e., it contains loops. An example of such a topology is shown in Fig. 5. Medium voltage grids often fall in an intermediate regime, with varying amounts of meshing.

Topology of a high voltage power grid. The visualization has been produced by Paul Cuffe using the methods described by Cuffe and Keane (2017) based on data of the NESTA archive. From Coffrin, Gordon, and Scott, 2014.

From a network science perspective a notable aspect of power grids is that they are geographically embedded. That is, nodes and lines have geographical locations and lengths. Line crossings are possible, but rare. These aspects are reflected in the connectivity or topology of the network, putting them in a different class than the ensembles typically studied in network science.

Further, as infrastructure has costs, power grids typically also contain as few lines as possible. As a consequence, power grids are typically sparse. The observed average degree for transmission grids falls between 2 and 5, and the degree distribution peaks between 2 and 3 and has an exponential tail (Solé et al., 2008; Wang, Scaglione, and Thomas, 2010).

B. Network ensembles and synthetic grid models

The statistical physics of networks often considers ensembles of networks to map out typical topological properties and to provide a reference class for actual real-world datasets. Power grids do not fully resemble any of these common network ensembles. In particular, they are not well captured by small-world networks (Cotilla-Sanchez et al., 2012). Although the clustering coefficient is similarly high, there is a much more pronounced connectivity at large scales. A Watts-Strogatz small-world network of such a low degree would likely end up being disconnected for moderately sized networks, and the nonlocal links in power grids tend to be much more clustered (Wang, Scaglione, and Thomas, 2010).

To capture the full topological structure of power grids, we need to define novel network ensembles: probability distributions on the space of networks that capture the notion of what looks like a power grid and what does not. This is typically done by specifying a random process that produces networks, either by rewiring links or by iteratively growing networks. Given such an embedded topology, typical power line parameters can then be used to derive the power grid’s properties (Table I).

Wang, Scaglione, and Thomas (2010) proposed a rewiring-based process with the aim of staying close in spirit to the construction of small-world networks while taking some specific identified properties of power grids into account. In contrast, growth-based models, starting with that of Schultz, Heitzig, and Kurths (2014b), grow the grid with random node placements and new connections. The model of Schultz, Heitzig, and Kurths (2014b) mimics the trade-off between global resilience and economy in the growth process and is able to recreate the exponential degree distribution. Stating these trade-offs explicitly allows one to study the trade-offs of dynamic and structural stability given by Plietzsch et al. (2016). A growth model for hierarchical networks was proposed and analyzed by Ódor and Hartmann (2018).

Soltan and Zussman (2016) considered node placement in addition to line generation. They considered average path length, clustering coefficient, the slope of node degree distributions, and the line length distributions and showed that their algorithm matched those in various North American power grids.

In these models the growth process makes use of topological and embedding information but does not consider the resulting energy flows in the growth stage. Birchfield et al. (2017) and Birchfield, Xu, and Overbye (2018) added consideration of power flows in the dc approximation and voltage profiles to the iterative growth of synthetic grids. Soltan, Loh, and Zussman (2019) fitted the growth process using a Gaussian mixture model. Espejo, Lumbreras, and Ramos (2019) presented a model that focuses instead on the historical plausibility of the growth process and economy versus robustness trade-offs.

Besides generating topologies from scratch, it is also possible to take an existing topology and vary it while keeping many aspects of it fixed. Such surrogate ensembles for spatial infrastructure networks were described by Wiedermann et al. (2016). A difficult open question is to find a fully satisfactory model for the cogeneration of a variety of networks (Hackl and Adey, 2019), or hierarchies of networks (Schultz et al., 2016).

C. Network datasets

The availability of real grid data is limited. Thus, a variety of test cases have been compiled for scientific use. Some are based on the available data of real-world grids. Some are synthetic ones trying to mimic the structural properties of real-world grids. In the following, we list resources that provide grid data without a claim of completeness.

• The classical IEEE test cases are probably the most heavily studied power grids of all, with some having been implemented decades ago. The datasets are partly synthetic and partly derived from real grids and can be obtained at the repository (Christie, 1999).

• Several online repositories have been created in recent years to provide larger, more diverse, and more recent test cases focusing on transmission grids (Birchfield, 2016; Pacific Northwest National Laboratory and National Rural Electric Cooperative Association, 2017; Farid et al., 2019) or distribution (“feeder”) grids (Schneider et al., 1991; Kavasseri and Ababei, 2021). In addition, several papers reviewed such repositories and the methodologies; see Coffrin, Gordon, and Scott (2014), Birchfield et al. (2017), and Schneider et al. (2018). These repositories again include both synthetic grids and approximate models of real grids, with a focus on North America. Large test grids with properties typical for European grids were provided by the Pegase project (Villella et al., 2012).

• The IEEE and Pegase test grids are included with others in the latest version of the software package matpower (Zimmerman, Murillo-Sanchez, and Thomas, 2011).

• Many test cases are restricted to static properties and do not include the parametrization necessary for dynamic simulations. To analyze the dynamics of these systems, some heuristic assumptions about present or future parameters are necessary. Realistic dynamic test cases mostly come in the form of small synthetic test cases developed by IEEE that can be found in repositories (Pacific Northwest National Laboratory and National Rural Electric Cooperative Association, 2017; Farid et al., 2019). ENTSO-E provides semisynthetic grid models mimicking the European grid (Semerow et al., 2015), but access is limited. A notable exception is the synthetic model of the northeastern U.S. published by Birchfield (2016), which includes a 25 000-node dynamic system, validated in the sense of the system of Xu, Birchfield, and Overbye (2018). A heuristic dynamical parametrization of the ENTSO-E–based system (Wiegmans, 2016) was recently given by Pagnier and Jacquod (2019c).

• Recently several initiatives started to map out power grids from publicly available data. A model of the German power grid was extracted from OpenStreetMaps (Medjroubi, Matke, and Kleinhans, 2015), and a model of the European grid was extracted from the ENTSO-E interactive grid map (Wiegmans, 2016). A variety of related datasets can be found on community websites (Open Energy Modelling Initiative, 2017; Open Power Systems Data, 2017). A collection of grid frequency time series was presented by Gorjão, Jumar et al. (2020).

V. Dynamics of Elementary Networks and Building Blocks

We introduce basic aspects of power system dynamics and stability for an elementary grid containing one transmission line. We first consider its static operation and its limitations before we turn to the dynamic stability.

A. Static solutions and voltage stability

We consider the steady state of an elementary circuit with one load node and one generator connected by a single transmission line to understand which factors limit the power transmission in different types of grids. The generator $s$ $s$ is treated as a slack node, so $V s = 1$ $V_{s} = 1$ and $δ s = 0$ $δ_{s} = 0$ are fixed; cf. Sec. III.A.3. Most load nodes draw the reactive power $Q n$ $Q_{n}$ at a fixed ratio to the real power, which is commonly specified in terms of the power factor $cos (ϑ)$ $\cos (ϑ)$ defined via $P n = - \sqrt P 2 n + Q 2 n cos (ϑ)$ $P_{n} = - \sqrt{P_{n}^{2} + Q_{n}^{2}} \cos (ϑ)$ and $Q n = - \sqrt P 2 n + Q 2 n sin (ϑ)$ $Q_{n} = - \sqrt{P_{n}^{2} + Q_{n}^{2}} \sin (ϑ)$ . We analyze the operation here as a function of the real power demand $P n < 0$ $P_{n} < 0$ .

To obtain the voltage magnitude $V n$ $V_{n}$ and phase angle $δ n$ $δ_{n}$ at the load node, we then solve the ac load flow equations (25), which are rewritten using the addition theorems of the sine and cosine

P n = G n n V 2 n + V n V s \sqrt B 2 n s + G 2 n s sin (δ n - γ n s), Q n = - B n n V 2 n - V n V s \sqrt B 2 n s + G 2 n s cos (δ n - γ n s),

$P_{n} = G_{n n} V_{n}^{2} + V_{n} V_{s} \sqrt{B_{n s}^{2} + G_{n s}^{2}} \sin (δ_{n} - γ_{n s}), Q_{n} = - B_{n n} V_{n}^{2} - V_{n} V_{s} \sqrt{B_{n s}^{2} + G_{n s}^{2}} \cos (δ_{n} - γ_{n s}),$

(91)

where the angle $γ n s$ $γ_{n s}$ is defined via

cos (γ n s) = + B n s \sqrt B 2 n s + G 2 n s, sin (γ n s) = - G n s \sqrt B 2 n s + G 2 n s .

$\cos (γ_{n s}) = \frac{+ B_{n s}}{\sqrt{B_{n s}^{2} + G_{n s}^{2}}}, \sin (γ_{n s}) = \frac{- G_{n s}}{\sqrt{B_{n s}^{2} + G_{n s}^{2}}} .$

We can then separate the equation for the voltage magnitude using $sin 2 + cos 2 = 1$ $\sin^{2} + \cos^{2} = 1$ , which yields

0 = (B 2 n n + G 2 n n) V 4 n + (P 2 n + Q 2 n) - (B 2 n s V 2 s + G 2 n s V 2 s - 2 B n n Q n + 2 G n n P n) V 2 n,

$0 = (B_{n n}^{2} + G_{n n}^{2}) V_{n}^{4} + (P_{n}^{2} + Q_{n}^{2}) - (B_{n s}^{2} V_{s}^{2} + G_{n s}^{2} V_{s}^{2} - 2 B_{n n} Q_{n} + 2 G_{n n} P_{n}) V_{n}^{2},$

(92)

which is easily solved for $V 2 n$ $V_{n}^{2}$ . The phase angle can then be obtained by solving one of Eqs. (91).

To understand the behavior of the system and the limitations to power transmission in more detail, we plot the voltage $V n$ $V_{n}$ as a function of the real power demand $| P n |$ $| P_{n} |$ in Fig. 6(a) for different values of the power factor. This plot is often referred to as the “nose curve” in power engineering (Machowski, Bialek, and Bumby, 2008). First, one observes that the quadratic equation (92) can have two solutions, where only the branch with the higher magnitude is relevant, as the other branch is unstable. Second, power transmission is generally limited. Physical solutions of Eq. (92) exist only if the real-power demand does not exceed the limit,

P 2 n \leq 14 B 2 n s V 4 s + B n s V 2 s Q n,

$P_{n}^{2} \leq \frac{1}{4} B_{n s}^{2} V_{s}^{4} + B_{n s} V_{s}^{2} Q_{n},$

(93)

which depends on the reactive power. The limit as well as the voltage magnitude $V n$ $V_{n}$ is lower when the load draws reactive power [ $Q n < 0$ $Q_{n} < 0$ , $tan (ϑ) > 0$ $\tan (ϑ) > 0$ ]. However, if the load node supplies reactive power instead [ $Q n > 0$ $Q_{n} > 0$ , $tan (ϑ) < 0$ $\tan (ϑ) < 0$ ], the voltage at the load node can be even larger than $V s$ $V_{s}$ and the transmittable real power is strongly enhanced. Grid operators hence aim to keep $tan (ϑ)$ $\tan (ϑ)$ small or even negative by supplying reactive power near the loads. This procedure is referred to as “reactive power compensation.”

Nose curve. (a) Voltage at the load node as a function of transmitted power for three different power factors: $tan (ϑ) = + 0.3, 0, - 0.3$ $\tan (ϑ) = + 0.3, 0, - 0.3$ (from left to right). (b) Voltage at the load node for $tan (ϑ) = 0$ $\tan (ϑ) = 0$ in the initial state (right) and after half of the connectivity is lost, i.e., $B n s$ $B_{n s}$ is halved (left). Parameters are $S base = 100 MW$ $S_{base} = 100 MW$ , $V s = 1 pu$ $V_{s} = 1 pu$ , and $B n s = 10 pu$ $B_{n s} = 10 pu$ . Solid (dashed) lines indicate stable (unstable) solutions of the load flow equations.

We further consider how a “voltage instability” can arise in elementary power grids. We assume that two parallel transmission lines exist, and that one of them is lost in a contingency situation. As a consequence, the effective susceptance $B n s$ $B_{n s}$ is halved and the nose curve contracts, as shown in Fig. 6(b). Further effects depend on the grid loading. If the real power demand $| P n |$ $| P_{n} |$ is not too high (such as the 200 MW value in the figure), a solution of the power flow equations still exists, albeit at a lower voltage. The consumers will thus experience a rapid drop in the voltage level. However, if the real power demand is too large (such as the 300 MW value in the figure), no static equilibrium solution exists anymore and a “voltage collapse” emerges in which voltages decline as governed by the power system dynamics. A comprehensive analysis of voltage stability in power grids was given by Van Cutsem and Vournas (2007).

B. Flow and voltage limits

Power transmission in current highest-voltage grids is typically limited by factors other than the ones discussed earlier. Most grids are constantly monitored and thoroughly regulated. Security limits have been formulated for current and voltage [see European Network of Transmission System Operators for Electricity (2004)], and emergency shutdowns are carried out when these limits are violated. Hence, the physical limits of power transmission are rarely met during normal operation.

In particular, extreme currents can lead to an Ohmic heating and eventually to a dangerous bending of overhead transmission lines due to thermal elongation. During the 2003 power outages in North America (U.S.-Canada Power System Outage Task Force, 2014) and Italy (Union for the Co-ordination of Transmission of Electricity, 2004) transmission lines hit trees, leading to a short-circuit fault. Hence, grid operators typically impose “thermal limits” for the currents and overloaded lines switch off automatically. Similarly, upper and lower limits for the voltage magnitude are imposed to guarantee power quality and avoid the danger of a complete voltage sack.

To see how these regulations limit power transmission in highest-voltage grids, we solve the load flow equations (91) as a function of $P n$ $P_{n}$ , assuming a fixed power factor set to a typical value of $cos (ϑ) = 0.95$ $\cos (ϑ) = 0.95$ , with $ϑ > 0$ $ϑ > 0$ . We use typical parameters for a 380 kV transmission line in Western Europe with $r / l = 0.03 Ω / km$ $r / l = 0.03 Ω / km$ and $x / l = 0.246 Ω / km$ $x / l = 0.246 Ω / km$ and a thermal limiting current of $I th = 2.58 kA$ $I_{th} = 2.58 kA$ (Oeding and Oswald, 2016), with $l$ $l$ the length of the line. We use the notation of Sec. V.A but give all quantities in physical units for better accessibility.

Figure 7 shows the results for two characteristic cases. We first consider a transmission line of length $l = 100 km$ $l = 100 km$ connecting the load node and a generator node, which has a fixed voltage magnitude $V g = 400 kV$ $V_{g} = 400 kV$ , i.e., slightly above the nominal voltage level of 380 kV. The current

| I n g | = ∣ ∣ ∣ | V n | e i δ n - | V g | r + i x ∣ ∣ ∣

$| I_{n g} | = | \frac{| V_{n} | e^{i δ_{n}} - | V_{g} |}{r + i x} |$

(94)

increases approximately linearly with the transmitted power. For $| P n | \geq 1.61 GW$ $| P_{n} | \geq 1.61 GW$ , the current exceeds the thermal line limit, which would lead to emergency measures, whereas the voltage magnitude stays well above the lower limit of $V lim = 0.9 \times 380 kV$ $V_{\lim} = 0.9 \times 380 kV$ . This illustrates that in densely populated areas with short lines, it is mainly the thermal limit for the current that limits the operation of a transmission line. In fact, we find a good estimate for the maximum transmittable real power using

P max \approx 3 380 kV \sqrt 3 \times I th = 1.698 GW .

$P_{\max} \approx 3 \frac{380 kV}{\sqrt{3}} \times I_{th} = 1.698 GW .$

(95)

Voltage stability can be an issue if generators are located far from loads. Figure 7 shows the current and voltage at the load node for a longer line with $l = 200 km$ $l = 200 km$ . The voltage limit $V lim$ $V_{\lim}$ is already hit for $| P n | \geq 1.01 GW$ $| P_{n} | \geq 1.01 GW$ , where the current is still well below the thermal limit.

Factors limiting real power transmission in ac power grids for two different characteristic cases (blue, short line with $l = 100 km$ $l = 100 km$ ; dashed red, long line with $l = 200 km$ $l = 200 km$ ). We consider a single transmission line connecting one power plant to a substation. Shown are the current $| I |$ $| I |$ , the voltage magnitude $| V n |$ $| V_{n} |$ , and the voltage angle $δ n$ $δ_{n}$ at the load node as a function of the real-power load $| P n |$ $| P_{n} |$ . The dotted line shows the limits $I th$ $I_{th}$ and $V lim$ $V_{\lim}$ for the current and the voltage.

Voltage limit violations can be handled to some extent by using reactive power compensation, as discussed in Sec. V.A. Decreasing the power angle $ϑ$ $ϑ$ at the load node allows for a higher power transmission. In contrast, the current limit represents a more severe limitation. A simple way to extend power transmission is dynamic line rating, where the thermal limiting current $I th$ $I_{th}$ is adapted to the ambient conditions (Douglass and Edris, 1996). On cold windy days, a higher current is acceptable, without the risk of overheating.

C. Dynamic stability

1. Local stability of a single generator

We now turn to the question of dynamical stability, starting with the most elementary setup: a single generator coupled to an infinite grid. The grid is assumed to be so large that it is not affected by generator dynamics, such that its voltage magnitude $V$ $V$ and phase angle $ϕ$ $ϕ$ are fixed. We consider here the question of angular stability, and thus assume that the generator voltage is constant at the steady-state value $E \circ$ $E^{\circ}$ . If the phenomenon of voltage collapse occurs, voltages drop dynamically such that a more refined treatment is necessary.

Assuming a lossless connection of the machine to the grid with the susceptance $x$ $x$ , the real power flow from the machine to the grid is given by $P el = K sin (δ - ϕ)$ $P_{el} = K \sin (δ - ϕ)$ and the swing equation (49) reads

2 H ω R ¨ δ + D ˙ δ = P mech - K sin (δ - ϕ),

$\frac{2 H}{ω_{R}} \ddot{δ} + D \dot{δ} = P_{mech} - K \sin (δ - ϕ),$

(96)

where $K = V E \circ / x$ $K = V E^{\circ} / x$ is an effective coupling strength. In normal operation the machine rotates with the frequency $ω R$ $ω_{R}$ and a fixed phase difference to the grid to supply a constant real power. Hence, we are interested in the fixed points of Eq. (96) given by

δ \circ + = ϕ + arcsin (P mech / K), δ \circ - = ϕ - arcsin (P mech / K) + π .

$δ_{+}^{\circ} = ϕ + \arcsin (P_{mech} / K), δ_{-}^{\circ} = ϕ - \arcsin (P_{mech} / K) + π .$

(97)

One finds that the $+$ $+$ solution is linearly stable, while the other branch is unstable. The solutions vanish in a saddle node on a circle bifurcation when $P mech = P crit = K$ $P_{mech} = P_{crit} = K$ . Hence, the power that the machine can receive and transmit to the grid is limited to a maximum value of $P crit$ $P_{crit}$ .

The nature of this bifurcation becomes more obvious with a mechanical analog. Consider a point particle with an effective mass $m eff = 2 H / ω R$ $m_{eff} = 2 H / ω_{R}$ and a friction coefficient $η = D$ $η = D$ moving in a tilted washboard potential

V eff (δ) = - P mech δ - K cos (δ - ϕ),

$V_{eff} (δ) = - P_{mech} δ - K \cos (δ - ϕ),$

(98)

where $δ$ $δ$ denotes the particle position. Newton’s equation for the particle $m eff ¨ δ = - η ˙ δ - \partial V eff / \partial δ$ $m_{eff} \ddot{δ} = - η \dot{δ} - \partial V_{eff} / \partial δ$ is then equivalent to the swing equation (96). The tilted washboard potential is shown in Fig. 8 for different values of tilting $P mech$ $P_{mech}$ . For low tilting, stable (unstable) fixed points exist in the local minima (maxima) of the potential (green and red dots). As the tilting increases, the minima and maxima approach each other and finally vanish for $P mech > P crit$ $P_{mech} > P_{crit}$ . The tilted washboard potential has been studied in great detail in statistical physics, particularly in the presence of noise; see Risken (1996) and references therein.

The tilted washboard potential [Eq. (98)] as a mechanical analog to the dynamics of a single generator in the second-order model for (a) $P mech = 0.1 K$ $P_{mech} = 0.1 K$ , (b) $P mech = 0.5 K$ $P_{mech} = 0.5 K$ , (c) $P mech = 0.9 K$ $P_{mech} = 0.9 K$ , (d) $P mech = 1.1 K$ $P_{mech} = 1.1 K$ . The potential $V eff$ $V_{eff}$ is given in units of K. Green disks, stable fixed points; red crosses, saddle points.

2. Bifurcations

To obtain a comprehensive understanding of the dynamics of a single generator, we also analyze its global stability properties. Given an initial state $δ (0), ˙ δ (0)$ $δ (0), \dot{δ} (0)$ , what happens for long times? Does the generator relax to the previously discussed stable fixed point?

A global phase-space portrait of the single machine system is shown in Fig. 9 for different values of $P mech / K$ $P_{mech} / K$ . One observes that in most cases a limit cycle exists that corresponds to a desynchronized generator. The generator cannot exchange real power with the grid in this state, as $P el = K sin (δ - ϕ)$ $P_{el} = K \sin (δ - ϕ)$ is oscillating and averages out to zero. The mechanical input power accelerates the generator until this input is balanced by the damping such that it has a higher frequency than the grid ( $˙ δ > 0$ $\dot{δ} > 0$ ).

Basin of attraction (red) of the stable fixed point of Eq. (96) in the $δ$ $δ$ - $˙ δ$ $\dot{δ}$ plane. Parameters are $2 H / ω R = 1$ $2 H / ω_{R} = 1$ , $D = 0.3$ $D = 0.3$ , $K = 8$ $K = 8$ , and (a) $P mech = 1$ $P_{mech} = 1$ , (b) $P mech = 2$ $P_{mech} = 2$ , (c) $P mech = 4$ $P_{mech} = 4$ , and (d) $P mech = 9$ $P_{mech} = 9$ . The white dots are the stable and unstable fix points. The black line gives the limit cycle.

We thus find three parameter regions with different global stability properties; see Figs. 9 and 10. For $P mech / K > 1$ $P_{mech} / K > 1$ no fixed point exists and the system will always converge to the limit cycle. If $P mech / K$ $P_{mech} / K$ is small or the damping $η$ $η$ is large, the system is globally stable: It converges to the attractive fixed point for almost all initial states. In the remaining part of the parameter space, the fixed points and the limit cycle coexist such that the long-term dynamics depends on the initial state.

Stability map of a single generator coupled to an infinite grid in the second-order model. Three different parameter regions exist: a globally stable fixed point, a globally attractive limit cycle, and a region of coexistence of a fixed point and a limit cycle. Bifurcations between these regions are shown as colored lines. The red dashed line is the approximation (102) for low friction. See Risken (1996) and Manik *et al.* (2014) for a more detailed analysis.

An analytical approximation for the border between the globally stable and the coexistence regime can be obtained in the low-friction limit using Lyapunov’s second method (Parks, 1992; Risken, 1996; Manik et al., 2014). We define an energy functional $E = m eff (˙ δ) 2 / 2 - K cos (δ)$ $E = m_{eff} (\dot{δ})^{2} / 2 - K \cos (δ)$ , setting $ϕ = 0$ $ϕ = 0$ for simplicity. Using the equation of motion for $δ$ $δ$ , we obtain

d E d t = m eff ¨ δ ˙ δ + K sin (δ) ˙ δ = P mech ˙ δ - η (˙ δ) 2 .

$\frac{d E}{d t} = m_{eff} \ddot{δ} \dot{δ} + K \sin (δ) \dot{δ} = P_{mech} \dot{δ} - η (\dot{δ})^{2} .$

(99)

If $E$ $E$ decreases on average for all initial conditions, the system is in the globally stable regime. The condition for the border between the globally stable and the coexistence regime is therefore obtained by setting $¯ ¯¯¯¯¯¯¯¯ ¯ d E / d t T = 0$ ${\bar{d E / d t}}^{T} = 0$ , where the bar denotes the average over one period $T$ $T$ . Evaluating this condition yields

P mech ¯ ˙ δ T = η ¯ ¯¯¯¯¯ ¯ (˙ δ) 2 T .

$P_{mech} {\bar{\dot{δ}}}^{T} = η {\bar{(\dot{δ})^{2}}}^{T} .$

(100)

We can now calculate $¯ ˙ δ T = 2 π / T$ ${\bar{\dot{δ}}}^{T} = 2 π / T$ and obtain

¯ ¯¯¯¯¯ ¯ (˙ δ) 2 T = 1 T \int T 0 (˙ δ) 2 d t = 1 T \int π - π ˙ δ d δ = 1 T \int π - π \sqrt 2 m eff \sqrt E (δ, ˙ δ) + K cos (δ) d δ .

${\bar{(\dot{δ})^{2}}}^{T} = \frac{1}{T} \int_{0}^{T} {(\dot{δ})}^{2} d t = \frac{1}{T} \int_{- π}^{π} \dot{δ} d δ = \frac{1}{T} \int_{- π}^{π} \sqrt{\frac{2}{m_{eff}}} \sqrt{E (δ, \dot{δ}) + K \cos (δ)} d δ .$

(101)

At the bifurcation, when a globally stable fixed point loses stability, there is a trajectory that satisfies $˙ δ = 0$ $\dot{δ} = 0$ at each successive peak of the potential landscape. This trajectory has $E peak = K$ $E_{peak} = K$ , which can be assumed to be constant in the low-friction limit. Replacing $E (δ, ˙ δ) \approx K$ $E (δ, \dot{δ}) \approx K$ , the integral in Eq. (101) can be evaluated in closed form and the border between the globally stable and the coexistence region [Eq. (100)] is given by

P mech = 4 η \sqrt K π \sqrt m eff .

$P_{mech} = \frac{4 η \sqrt{K}}{π \sqrt{m_{eff}}} .$

(102)

Equation (102) matches the results from a direct numerical evaluation for small values of $η$ $η$ , as shown in Fig. 10. For large values of $η$ $η$ it slightly overestimates the globally stable parameter region.

3. Probabilistic stability measures

Elementary systems allow for a comprehensive understanding of the bifurcation structure and global stability. However, this approach is no longer feasible for high-dimensional networked systems, such that more versatile methods are needed. Recently, probabilistic approaches have attracted strong interest (Menck et al., 2013). For instance, one can quantify the stability of a dynamical system by the probability of returning to the desired fixed point after a random perturbation. The set of all initial states guaranteeing convergence, the basin of attraction, is shown in Fig. 9 for a single generator and different system parameters. While the geometry of the basin can be rather intricate in higher dimensions, its size is easy to measure by randomly drawing a sufficient number of initial states from a predefined set $A$ $A$ and numerically simulating their dynamics. It is then the same as the following Monte Carlo integration of the indicator function $1 B$ $1_{B}$ on the basin of attraction:

^P = P (x (\infty) = x * | x (0) \in A) = 1 | A | \int A 1 B (x) d x .

$\hat{P} = P (x (\infty) = x^{*} | x (0) \in A) = \frac{1}{| A |} \int_{A} 1_{B} (x) d x .$

(103)

Different aspects of stability can be understood by choosing an appropriate set $A$ $A$ . For instance, restricting perturbations to a single node reveals weak spots in the grid (Menck et al., 2014).

The process of sampling such a probability is a Bernoulli process. Its sampling uncertainty decreases with the number of samples $N$ $N$ , which is asymptotically given by the binomial confidence interval

^P \pm \sqrt^P (1 -^P) / N,

$\hat{P} \pm \sqrt{\hat{P} (1 - \hat{P}) / N},$

(104)

although better estimates are possible for small values of $N$ $N$ (Agresti and Coull, 1998). Notably the number of samples needed does not increase with the dimension of the problem, whereas the individual samples are typically more expensive to obtain.

The asymptotic convergence to a fixed point is an important aspect of stability, but it is not sufficient in the context of power grids. One can generalize the current approach, requiring that a trajectory does not exceed certain operational limits or, equivalently, that it does not not leave a “survival” region $S$ $S$ . This leads to the following definition of the survivability (Hellmann et al., 2016):

P (x (t) \in S \forall t | x (0) \in A) .

$P (x (t) \in S \forall t | x (0) \in A) .$

(105)

Survivability requires the specification of a survival region $S$ $S$ , as shown in Fig. 11 for a single generator, assuming that the survival region is given by $| ˙ δ | \leq \pm 5 rad / s$ $| \dot{δ} | \leq \pm 5 rad / s$ . This assumption is useful in the study of power grids, as generators have many protection circuits that would switch them off to avoid mechanical damage. Probabilistic measures have also been generalized to cover the case of repeated perturbations and stochastic systems (Schultz et al., 2018; Lindner and Hellmann, 2019). Furthermore, the uniform sampling from the set $A$ $A$ can be replaced by more general distributions of initial states.

Survivability region (red) of the stable fixed point of Eq. (96), with the survival region given by the area between the lines $˙ δ = \pm 5 rad / s$ $\dot{δ} = \pm 5 rad / s$ . The basin of attraction of Fig. 9 is shaded in orange. The parameters are $2 H / ω R = 1$ $2 H / ω_{R} = 1$ , $D = 0.3$ $D = 0.3$ , $K = 8$ $K = 8$ , and (a) $P = 1$ $P = 1$ , (b) $P = 2$ $P = 2$ , (c) $P = 4$ $P = 4$ , and (d) $P = 9$ $P = 9$ . The white dots are the stable and unstable fixed points. The black line gives the limit cycle.

4. Stochastic stability

Thus far we have analyzed the stability of the swing equation for fixed system parameters. But real-world systems will typically face disturbances and noise. How does this affect the stability of the system?

We now analyze how the stability of the swing equation (96) is affected if the mechanical input power of a machine fluctuates in time as

P mech (t) = ¯ P + ξ (t) .

$P_{mech} (t) = \bar{P} + ξ (t) .$

(106)

Notably this problem can be solved almost fully analytically in the case of white noise $ξ (t)$ $ξ (t)$ using Kramers’s escape rate theory (Risken, 1996). Using the analogy to a particle moving in a tilted washboard potential [Eq. (98)], the main question is as follows: Can the particle overcome the potential barrier and thus escape the vicinity of the stable fixed point as illustrated in Figs. 12(a) and 12(b)? What is the probability of such a desynchronization event?

Desynchronization of the swing equation due to a fluctuating power input. (a),(b) When the input power $P mech$ $P_{mech}$ fluctuates, the generator can lose synchrony to the grid after an escape time $τ$ $τ$ . (c) Kramers’s escape rate theory predicts the escape process for Gaussian white noise. The theoretical prediction [Eq. (107) (black lines)] accurately predicts the mean escape times $¯ τ$ $\bar{τ}$ obtained from direct numerical simulations averaged over 2000 sample processes (symbols). Results are shown for an extremely loaded connection with $K = 1$ $K = 1$ and $P = 0.95$ $P = 0.95$ and a damping coefficient of $D = 8 \times 10 - 5 s - 1 \times ω R$ $D = 8 \times 10^{- 5} s^{- 1} \times ω_{R}$ (intermediate damping). From Schäfer *et al.*, 2017.

Kramers’s theory gives the average time $¯ τ$ $\bar{τ}$ until a particle escapes (van Kampen, 2007; Schäfer et al., 2017) as follows:

¯ τ = 2 π ϱ \sqrt V'' eff (δ min) | V'' eff (δ max) | exp (2 η Δ V eff σ 2),

$\bar{τ} = \frac{2 π ϱ}{\sqrt{V_{eff}^{''} (δ_{\min}) | V_{eff}^{''} (δ_{\max}) |}} \exp (\frac{2 η Δ V_{eff}}{σ^{2}}),$

(107)

where $σ 2$ $σ^{2}$ is the variance of the Gaussian white noise, $η = D$ $η = D$ is the effective friction, and $2 ϱ = η + \sqrt η 2 + (8 H / ω R) | V'' eff (δ max) |$ $2 ϱ = η + \sqrt{η^{2} + (8 H / ω_{R}) | V_{eff}^{''} (δ_{\max}) |}$ for intermediate damping. The quantities $V'' (δ min, max)$ $V^{''} (δ_{\min, \max})$ are the second derivatives of the potential evaluated at the local minimum and maximum $δ min, max$ $δ_{\min, \max}$ , respectively, and

Δ V eff = V eff (δ max) - V eff (δ min)

$Δ V_{eff} = V_{eff} (δ_{\max}) - V_{eff} (δ_{\min})$

(108)

is the potential barrier that the particle has to overcome. Numerical simulations show an excellent agreement with this prediction, as indicated in Fig. 12(c).

Kramers’s formula (107) reveals how the system parameters affect the stability. The main factors are noise strength $σ$ $σ$ , damping parameter $D$ $D$ , and potential barrier $Δ V eff$ $Δ V_{eff}$ , which enter Kramers’s formula exponentially. In contrast, $¯ τ$ $\bar{τ}$ depends less sensitively on the inertia $H$ $H$ , which enters only algebraically.

The barrier is determined by the effective tilting of the potential, i.e., the loading of the grid, and vanishes as $¯ P \to P crit$ $\bar{P} \to P_{crit}$ . Figure 12 demonstrates the escape process for an extreme load, whereas $¯ τ$ $\bar{τ}$ is longer by orders of magnitude for realistic loading levels. However, $¯ τ$ $\bar{τ}$ can be drastically reduced if the noise is no longer Gaussian but instead intermittent, which is typical for wind power systems (Schmietendorf, Peinke, and Kamps, 2017). In larger power grids, different escape routes can exist, corresponding to different parts of the grid loosing synchrony (Schäfer et al., 2017). In Sec. VIII.D we reconsider the stability of the swing equation under several generalizations.

VI. Synchronization and Steady States in Complex Grids

A. The need for synchrony

The stable steady operation of a power grid requires perfect phase locking of all synchronous machines throughout the grid (Rohden et al., 2012; Dörfler, Chertkov, and Bullo, 2013; Motter et al., 2013). That is, all machines in a connected network have to run at the same frequency with a well-defined phase difference between them; otherwise, no steady power flow is possible. From Eq. (73) we find that if the nodes operate at their nominal voltage amplitude $E j = 1$ $E_{j} = 1$ , the real power flow between two nodes $i$ $i$ and $j$ $j$ is determined by their phase difference

P i \to j = K i j {sin (γ i j) + sin [δ i (t) - δ j (t) - γ i j]},

$P_{i \to j} = K_{i j} {\sin (γ_{i j}) + \sin [δ_{i} (t) - δ_{j} (t) - γ_{i j}]},$

(109)

where $Y eff i j = K i j e i (γ i j + π / 2)$ $Y_{i j}^{eff} = K_{i j} e^{i (γ_{i j} + π / 2)}$ is the effective nodal admittance and the term $sin (γ i j)$ $\sin (γ_{i j})$ comes from the diagonal of $Y eff$ $Y^{eff}$ . The same expression is used for the aggregated model introduced in Sec. III.C.9 with $γ i j = 0$ $γ_{i j} = 0$ .

Under fault conditions, a fast transient instability can occur when a group of machines accelerates relative to the others. The relative phases $δ i, j (t)$ $δ_{i, j} (t)$ between the groups are then no longer locked and there is no steady power flow $P i \to j$ $P_{i \to j}$ . If the fault clears, another transient may either diverge or converge. Transient instabilities on timescales of cycles typically lead to the shutdown (“tripping”) of both transmission lines and generators.

Furthermore, violations of perfect phase locking are observed during “interarea oscillations,” where electric power oscillates across the grid at low frequencies of 0.1–10 Hz (Klein, Rogers, and Kundur, 1991; Rogers, 2012). Inter-area oscillations repeatedly occur after disturbances such as the loss of a generator and are typically damped out in minutes. However, oscillations may also grow in exceptional cases, eventually leading to a system split or blackout. For instance, a cascade of failures led to the tripping of several transmission lines and the loss of the stable steady state during the 2006 European power outage (Union for the Coordination of Transmission of Electricity, 2007). Oscillations grew rapidly and the grid finally separated into three mutually asynchronous areas.

We are thus led to a fundamental question of power system stability: Do the equations of motion admit a stable solution where all nodes are perfectly phase locked,

δ j (t) = Ω t + δ \circ j for j \in {1, \dots, N},

$δ_{j} (t) = Ω t + δ_{j}^{\circ} for j \in {1, \dots, N},$

(110)

given the network topology and the power injections? Synchrony is about phase dynamics; therefore, it is natural to focus on the phases as dynamic variables in the following. To simplify the analysis, we neglect voltage dynamics for the time being and assume that the phase dynamics at the nodes is well described by the swing equation (49). As previously discussed, this equation is generic and arises similarly for single synchronous machines, grid-forming inverters with power control, and aggregated dynamical models. Absorbing the various factors in front of $¨ δ i$ ${\ddot{δ}}_{i}$ and $˙ δ i$ ${\dot{δ}}_{i}$ into reduced inertia $J i$ $J_{i}$ and damping $D i$ $D_{i}$ coefficients, the central equation in question then reads

J i ¨ δ i = P i - D i ˙ δ i - \sum j K i j [sin (γ i j) + sin (δ i - δ j - γ i j)] .

$J_{i} {\ddot{δ}}_{i} = P_{i} - D_{i} {\dot{δ}}_{i} - \sum_{j} K_{i j} [\sin (γ_{i j}) + \sin (δ_{i} - δ_{j} - γ_{i j})] .$

(111)

The general problem of synchronization and dynamical stability with voltage dynamics is more involved, and we present an outlook at the end of this section.

Versions of this model have also been studied in the theoretical physics literature under the name of the second-order Kuramoto model; cf. Chap. 7 of Rodrigues et al. (2016). However, note that different aspects of synchronization have traditionally attracted more attention in theoretical physics. Starting from the seminal work of Kuramoto (1975) many researchers have investigated the transition from incoherence to partial synchronization, where some nodes lock their frequencies. This is not sufficient for power systems. Furthermore, many essential results on partial synchronization have been obtained in a large- $N$ $N$ mean-field limit. A generalized large- $N$ $N$ limit for complex networks and the problem of full synchronization were treated only recently by Kuehn and Throm (2019).

B. Linear stability of the phase dynamics

We begin by assuming that a solution to the static equations is given. We then study whether its phase dynamics is stable to small perturbations. Consider a small perturbation around the synchronous fixed point $δ i (t) = Ω t + δ \circ i + α i (t)$ $δ_{i} (t) = Ω t + δ_{i}^{\circ} + α_{i} (t)$ . To linear order, perturbations evolve according to

J i ¨ α i = - D i ˙ α i - N \sum j = 1 Λ i j α j,

$J_{i} {\ddot{α}}_{i} = - D_{i} {\dot{α}}_{i} - \sum_{j = 1}^{N} Λ_{i j} α_{j},$

(112)

with the matrix $Λ \in R N \times N$ $Λ \in R^{N \times N}$ defined as

Λ i j = ⎧ ⎪ ⎨ ⎪ ⎩ - K i j cos (δ \circ i - δ \circ j - γ i j) if i \neq j, \sum n \neq i K i n cos (δ \circ i - δ \circ n - γ i j) if i = j .

$Λ_{i j} = {\begin{cases} - K_{i j} \cos (δ_{i}^{\circ} - δ_{j}^{\circ} - γ_{i j}) & if i \neq j, \\ \sum_{n \neq i} K_{i n} \cos (δ_{i}^{\circ} - δ_{n}^{\circ} - γ_{i j}) & if i = j . \end{cases}$

(113)

The fixed point is linearly stable if all perturbations $α$ $α$ are damped. We note that power grid synchrony is not affected if all angles $δ i$ $δ_{i}$ are shifted by a constant. Thus, this mode is excluded from the stability analysis. One can show that stability is primarily determined by $Λ$ $Λ$ : The synchronous state is linearly stable if all eigenvalues of $Λ$ $Λ$ are positive except for the one corresponding to a global shift of all phases. This criterion becomes especially simple if line losses can be neglected. $γ i j = 0$ $γ_{i j} = 0$ , and $Λ$ $Λ$ is a weighted, symmetric Laplacian matrix; cf. Sec. II. Stability is guaranteed as long as the “weights” $K i j cos (δ \circ i - δ \circ j)$ $K_{i j} \cos (δ_{i}^{\circ} - δ_{j}^{\circ})$ are positive for all links $(i, j)$ $(i, j)$ ; i.e., if no line is overloaded,

| δ \circ i - δ \circ j | < π / 2 .

$| δ_{i}^{\circ} - δ_{j}^{\circ} | < π / 2 .$

(114)

This criterion is sufficient, but generally not necessary. It can be viewed as a generalization of the stability criterion for a single generator discussed in Sec. V.C.1.

The situation becomes more involved when line losses cannot be neglected. The condition $| δ \circ i - δ \circ j - γ i j | < π / 2$ $| δ_{i}^{\circ} - δ_{j}^{\circ} - γ_{i j} | < π / 2$ is then not sufficient and more subtle necessary conditions are derived (Skar, 1980). Further effective criteria for the definiteness of signed Laplacians were given by Song, Hill, and Liu (2015) and Chen et al. (2016). Lyapunov exponents and vectors for this system were studied by Bosetti and Khan (2018). Note that phase differences in transmission grids are typically much smaller. Other factors limit power transmission, at least for short transmission lines; see Fig. 7.

Stability can be lost if a parameter of the network is varied, e.g., if power injections and grid loads increase or transmission lines fail (Coletta and Jacquod, 2020). Eventually, the second Laplacian eigenvalue tends to zero and the fixed point is lost in an inverse saddle-node bifurcation on a circle (Manik et al., 2014). Notably bifurcation sets and linear stability properties become much more intriguing if voltage dynamics is included (Schmietendorf et al., 2014; Ma et al., 2016; Sharafutdinov et al., 2018) or in inverter-based grids (Schiffer et al., 2014). An overview of different aspects of dynamical stability was provided by Gajduk, Todorovski, and Kocarev (2014).

C. Existence, multiplicity, and properties of synchronous states

We now return to the question of whether a stable synchronous state exists for a given power grid. We first review some important properties of synchronous states and then review a theoretical approach to systematically compute and classify all such states.

1. Properties of synchronous states

The proper operation of a power grid requires perfect phase locking according to Eq. (110). Substituting into the equations of motions (111) yields the following fundamental condition for the phases and active power flows in the synchronous state of a power grid:

0 = P i - D i Ω - \sum j K i j [sin (γ i j) + sin (δ \circ i - δ \circ j - γ i j)] .

$0 = P_{i} - D_{i} Ω - \sum_{j} K_{i j} [\sin (γ_{i j}) + \sin (δ_{i}^{\circ} - δ_{j}^{\circ} - γ_{i j})] .$

(115)

We note that the flows are asymmetric for $γ \neq 0$ $γ \neq 0$ due to Ohmic losses. In fact, one can express the losses at a transmission line $(i, j)$ $(i, j)$ as

P loss i j = P i \to j + P j \to i = 2 K i j sin (γ i j) [1 - cos (δ \circ i - δ \circ j)] .

$P_{i j}^{loss} = P_{i \to j} + P_{j \to i} = 2 K_{i j} \sin (γ_{i j}) [1 - \cos (δ_{i}^{\circ} - δ_{j}^{\circ})] .$

Losses vanish when $γ i j = 0$ $γ_{i j} = 0$ and when there is no power flowing on the line, i.e., when $δ \circ i = δ \circ j$ $δ_{i}^{\circ} = δ_{j}^{\circ}$ .

Synchrony does not necessarily imply that the grid operates exactly at the reference frequency of 50 or $60 Hz$ $60 Hz$ . In fact, we can determine the equilibrium frequency directly from Eq. (115). Summing over all nodes $i \in {1, \dots, N}$ $i \in {1, \dots, N}$ and solving for $Ω$ $Ω$ yield

Ω = \sum i P i - \sum i < j P loss i j \sum i D i .

$Ω = \frac{\sum_{i} P_{i} - \sum_{i < j} P_{i j}^{loss}}{\sum_{i} D_{i}} .$

(116)

We recall that the effective parameters $D i$ $D_{i}$ include both damping and the action of primary load-frequency control, such that the current result is equivalent to Eq. (89). The deviation from the reference frequency $Ω$ $Ω$ is given by the ratio of the power imbalance, including losses, and the cumulative primary control strength. We return to the analysis of the bulk frequency in Sec. VIII.A.

2. Phase cohesiveness and necessary conditions

We are particularly interested in synchronous states where angle differences are small; cf. the discussion in Sec. V.A. Small angle differences are also essential for many analytic results on the stability of the swing equation; cf. the linear stability condition (114). To state this aspect more precisely we define the cohesiveness as follows: A state is called phase cohesive if all phases $δ \circ i$ $δ_{i}^{\circ}$ lie in an arc of length $ζ \in [0, π / 2)$ $ζ \in [0, π / 2)$ denoted as $A (ζ)$ $A (ζ)$ .

The definition can be applied to deriving necessary conditions for synchronization. For simplicity, assume that the power in the grid is balanced and that we have no losses $\sum i P i = 0$ $\sum_{i} P_{i} = 0$ , $γ i j = 0$ $γ_{i j} = 0$ . Hence, we have $Ω = 0$ $Ω = 0$ in Eq. (110) and the condition for phase locking in the second-order power grid model (111) reads

P i = N \sum n = 1 K i n sin (δ \circ i - δ \circ n) .

$P_{i} = \sum_{n = 1}^{N} K_{i n} \sin (δ_{i}^{\circ} - δ_{n}^{\circ}) .$

(117)

For a phase cohesive sync state, the magnitude of the sum can be bounded from above by $sin (ζ) \sum n K i n$ $\sin (ζ) \sum_{n} K_{i n}$ such that Eq. (117) can be satisfied only if, for all $i \in {1, \dots, N}$ $i \in {1, \dots, N}$ ,

| P i | \leq sin (ζ) N \sum n = 1 K i n .

$| P_{i} | \leq \sin (ζ) \sum_{n = 1}^{N} K_{i n} .$

(118)

Furthermore, we can consider any two nodes $i$ $i$ and $j$ $j$ of the grid and evaluate the difference of the conditions (117). We can again bound the sums on the right-hand side which yields the following condition for all $i, j \in {1, \dots, N}$ $i, j \in {1, \dots, N}$ :

| P i - P j | \leq sin (ζ) N \sum n = 1 K i n + K j n .

$| P_{i} - P_{j} | \leq \sin (ζ) \sum_{n = 1}^{N} K_{i n} + K_{j n} .$

(119)

The two conditions are necessary for a cohesive sync state to exist, but they are not sufficient by any means. Nevertheless, they reveal two important facts: (i) Synchronization requires a sufficient connectivity of the grid, while a strong divergence of the $P i$ $P_{i}$ impedes synchronization. (ii) If the phase divergence $ζ$ $ζ$ is reduced, we have to increase the coupling or decrease the divergence of $P i$ $P_{i}$ . Further generalizations and refinements were discussed by Chopra and Spong (2009), Ainsworth and Grijalva (2013), and Dörfler, Chertkov, and Bullo (2013).

3. Existence of solutions and multistability

We now turn to the question of whether solutions to the power flow equations exist. Notably power grids can admit multiple stable and unstable states of operation, even when power injections are fixed. Sudden changes of the power injections or the grid topology can trigger transitions between these states, with drastic consequences for grid operation and stability. Circulating power flows can emerge, which are generally undesired, as they increase line loads and Ohmic losses. A classic example of such circulating power flows is the Lake Erie loop in North America (Coletta et al., 2016; Jafarpour et al., 2019). If the system finds itself in a situation in which no stable synchronous state exists, system collapse is an inevitable consequence. We begin with the simplest case of a power grid without Ohmic losses such that the fixed points are determined by Eq. (117). The central questions now are as follows: When does this equation have multiple solutions that correspond to a dynamically stable power flow, and why?

We start with an elementary example: a simple ring network where $P i = 0$ $P_{i} = 0$ for all nodes $i \in {1, \dots, N}$ $i \in {1, \dots, N}$ and homogeneous line parameters $K$ $K$ (Schröder, Timme, and Witthaut, 2017). One can easily verify that the fixed points of this network are given by $δ * i = 2 π i m / N$ $δ_{i}^{*} = 2 π i m / N$ , with a parameter $m \in {- N / 2, - N / 2 + 1, \dots, + N / 2}$ $m \in {- N / 2, - N / 2 + 1, \dots, + N / 2}$ . Fixed points with $m \in (- N / 4, + N / 4)$ $m \in (- N / 4, + N / 4)$ satisfy Eq. (114) and are guaranteed to be dynamically stable; all other fixed points are potentially unstable. This simple example already suggests the fundamental features of multiplicity of sync states. Power flow solutions in lossless grids can differ in cycle flows. In fact, each fixed point corresponds to a flow $f = K sin (2 π m / N)$ $f = K \sin (2 π m / N)$ around the cycle. The number of stable fixed points increases with the size of the loop and is limited to one if $N \leq 4$ $N \leq 4$ .

Using these insights, we introduce a general method to construct all stable fixed points of a generic lossless power grid following Manik, Timme, and Witthaut (2017). The main idea is to shift the focus from nodes to the edges and cycles of the network. Hence, we define a vector of flows on the network’s edges as

F = B d sin (E ⊤ δ) \in R L,

$F = B_{d} \sin (E^{⊤} δ) \in R^{L},$

(120)

where $E$ $E$ denotes the node-edge incidence matrix (4) and the sine function is taken elementwise. The steady-state conditions (117) then become

P = E F .

$P = E F .$

(121)

Fixed points are now obtained via a stepwise procedure: We first construct all solutions of Eq. (121) and then reject those that are incompatible with Eq. (120).

The first step is based on the following observation: The kernel of the incidence matrix $E$ $E$ corresponds exactly to cycle flows in the network. Hence, the general solution of Eq. (121) can be written as

F = F (s) + C f,

$F = F^{(s)} + C f,$

(122)

where $F (s) \in R L$ $F^{(s)} \in R^{L}$ is a specific solution, $f \in R L - N + 1$ $f \in R^{L - N + 1}$ gives the strength of the cycle flows, and the matrix $C$ $C$ is the cycle-edge incidence matrix (5). For each flow vector $F$ $F$ , we can construct the associated nodal phases as follows. Start with the slack node $s$ $s$ and set $δ \circ s = 0$ $δ_{s}^{\circ} = 0$ . Proceed to a neighboring node $j$ $j$ . If one assumes that the connecting edge $e^= (j, s)$ $e \hat{=} (j, s)$ is oriented from node $s$ $s$ to node $j$ $j$ , the phase value reads $δ \circ j = δ \circ s + Δ e$ $δ_{j}^{\circ} = δ_{s}^{\circ} + Δ_{e}$ , where the phase difference $Δ e$ $Δ_{e}$ is reconstructed from the flow $F e$ $F_{e}$ by inverting Eq. (120), which yields the two following possible solutions:

Δ + e = arcsin (F e / K e), Δ - e = π - arcsin (F e / K e) .

$Δ_{e}^{+} = \arcsin (F_{e} / K_{e}), Δ_{e}^{-} = π - \arcsin (F_{e} / K_{e}) .$

(123)

We must choose one of the solutions, and we keep track of this choice by decomposing the edge set as $E = E + \cup E -$ $E = E_{+} \cup E_{-}$ , where $E \pm = {e \in E | Δ e = Δ \pm e}$ $E_{\pm} = {e \in E | Δ_{e} = Δ_{e}^{\pm}}$ . Not all solutions obtained in this way are physically correct. We can get the physically correct ones by making sure that the sum of the phase differences around any fundamental cycle yields 0 or an integer multiple of $2 π$ $2 π$ , which is expressed via the following winding number condition:

ϖ c = 1 2 π L \sum e = 1 C e, c Δ \pm e \in Z .

$ϖ_{c} = \frac{1}{2 π} \sum_{e = 1}^{L} C_{e, c} Δ_{e}^{\pm} \in Z .$

(124)

Fixed points with $E - = \emptyset$ $E_{-} = \emptyset$ satisfy Eq. (114) and are always stable, while states with $E - \neq \emptyset$ $E_{-} \neq \emptyset$ are typically unstable (Manik et al., 2014; Delabays, Coletta, and Jacquod, 2017; Manik, Timme, and Witthaut, 2017).

The cycle flow approach yields numerous analytical insights into the occurrence of multiplicity of sync states:

• No multiplicity in trees.—In a lossless tree network, either there is no fixed point or there are $2 N - 1$ $2^{N - 1}$ fixed points, of which one is stable and $2 N - 1$ $2^{N - 1}$ are unstable (Manik, Timme, and Witthaut, 2017).

• Similarly, multiplicity is ruled out in dense networks, as the fundamental cycles are too small (Taylor, 2012). Hence, all fixed points except one include edges with $| δ \circ i - δ \circ j | > π / 2$ $| δ_{i}^{\circ} - δ_{j}^{\circ} | > π / 2$ .

• In a simple ring network, one can derive explicit upper and lower bounds for the number of fixed points. Generally, the number of fixed points increases with the size of the ring and decreases with the loading of the grid (Ochab and Gora, 2010; Delabays, Coletta, and Jacquod, 2016; Manik, Timme, and Witthaut, 2017).

• In plane networks, one can show that the winding numbers are unique. That is, two fixed points cannot have the same winding numbers and the same decomposition $E = E + \cup E -$ $E = E_{+} \cup E_{-}$ . This can be used to derive upper bounds for the number of stable fixed points (Delabays, Coletta, and Jacquod, 2017; Manik, Timme, and Witthaut, 2017).

• One can derive estimates for the number of stable fixed points using scaling relations (Manik, Timme, and Witthaut, 2017); an example is shown in Fig. 13. Further numerical results for large sparse networks were discussed by Mehta et al. (2015) and Xi, Dubbeldam, and Lin (2017).

Multiplicity of sync states: scaling of the number of stable fixed points in a grid with two cycles. The grid consists of two rings with $n + 1$ $n + 1$ nodes sharing one edge such that $N = 2 n$ $N = 2 n$ . Edge weights are assumed to be homogeneous and loads vanish ( $P i = 0$ $P_{i} = 0$ ). The analytic scaling result $0.1576 (n 2 + 2 n)$ $0.1576 (n^{2} + 2 n)$ (solid line) matches the numerically exact results (circles) well. Adapted from Manik, Timme, and Witthaut, 2017.

The previous approach can be generalized to include Ohmic losses (Balestra et al., 2019). One again first focuses on the lines and considers the flows $F i j$ $F_{i j}$ and the losses $L i j = g j k [1 - cos (δ \circ i - δ \circ j)]$ $L_{i j} = g_{j k} [1 - \cos (δ_{i}^{\circ} - δ_{j}^{\circ})]$ as basic variables. The power balance equation is linear in these variables, and one can construct the general solution of this set of equations. Among this large set of solution candidates, the correct solutions are found by imposing the winding number conditions (124) and the additional constraint

(F i j b i j) 2 + (L i j g i j - 1) 2 = 1 for all edges (i, j),

${(\frac{F_{i j}}{b_{i j}})}^{2} + {(\frac{L_{i j}}{g_{i j}} - 1)}^{2} = 1 for all edges (i, j),$

(125)

which follows from $sin 2 + cos 2 = 1$ $\sin^{2} + \cos^{2} = 1$ . One then finds that Ohmic losses can have two conflicting effects on the existence and number of stable steady states. On the one hand, high losses must be compensated for by higher power flows, which may decrease their number. On the other hand, Ohmic losses can stabilize certain solution branches and thus foster multistability.

D. Nonlinear stability and explicit synchronization criteria

The existence of a phase-locked state is a prerequisite for the stable operation of a power grid. But even in elementary networks phase-locked states and limit cycles can coexist (see Fig. 9), and whether or not the phases lock depends on the initial conditions. In the following, we review explicit criteria for synchronization based on the dynamics (111) and discuss the main factors determining a network’s synchronization capability. A comprehensive overview is provided in the review by Dörfler and Bullo (2014).

1. Sufficient criteria for dense networks

We first present a sufficient condition for dense networks, which arise naturally in power system models after Kron reduction; cf. Sec. III.C.6. We are interested mainly in second-order power grid model dynamics [Eq. (111)], but it is simpler to start with the following overdamped limit:

D i ˙ δ i = P i + \sum j K i j [sin (γ i j) - sin (δ i - δ j - γ i j)] .

$D_{i} {\dot{δ}}_{i} = P_{i} + \sum_{j} K_{i j} [\sin (γ_{i j}) - \sin (δ_{i} - δ_{j} - γ_{i j})] .$

(126)

Dörfler and Bullo (2012b) proved that this system achieves phase cohesiveness and frequency synchronization if the coupling is strong compared to the differences in natural frequencies [ $ω nat i = P i + \sum j K i j sin (γ i j)$ $ω_{i}^{nat} = P_{i} + \sum_{j} K_{i j} \sin (γ_{i j})$ ]. They quantified these two opposing forces in terms of the parameters

Γ min ≔ N min i \neq j (K i j D i cos (γ i j)),

$Γ_{\min} ≔ N \min_{i \neq j} (\frac{K_{i j}}{D_{i}} \cos (γ_{i j})),$

(127)

Γ crit ≔ 1 cos (γ max) (max i \neq j ∣ ∣ ∣ ∣ ω nat i D i - ω nat j D j ∣ ∣ ∣ ∣ + 2 max i N \sum j = 1 K i j D i sin (γ i j)) .

$Γ_{crit} ≔ \frac{1}{\cos (γ_{\max})} (\max_{i \neq j} | \frac{ω_{i}^{nat}}{D_{i}} - \frac{ω_{j}^{nat}}{D_{j}} | + 2 \max_{i} \sum_{j = 1}^{N} \frac{K_{i j}}{D_{i}} \sin (γ_{i j})) .$

(128)

Furthermore, they defined $ζ min \in [0, π / 2 - γ max)$ $ζ_{\min} \in [0, π / 2 - γ_{\max})$ and $ζ max \in (π / 2, π]$ $ζ_{\max} \in (π / 2, π]$ as the solutions of the equation $sin (ζ min) = sin (ζ max) = cos (γ max) Γ crit / Γ min$ $\sin (ζ_{\min}) = \sin (ζ_{\max}) = \cos (γ_{\max}) Γ_{crit} / Γ_{\min}$ whenever this solution exists.

Then they derived the following theorem: If $Γ min > Γ crit$ $Γ_{\min} > Γ_{crit}$ and the angles are initially not too different [all $δ i (0)$ $δ_{i} (0)$ are in the arc $A (ζ max)$ $A (ζ_{\max})$ ], the network will achieve exponential frequency synchronization; i.e., all frequencies $˙ δ i$ ${\dot{δ}}_{i}$ converge exponentially fast to a common frequency $Ω$ $Ω$ . The phases remain cohesive and all $δ i (t)$ $δ_{i} (t)$ reach the arc $¯ A (ζ min)$ $\bar{A} (ζ_{\min})$ .

Dörfler and Bullo further generalized this result to the second-order power system model (111) if the damping was high enough. More precisely, if the ratio $ε = J max / D min$ $ε = J_{\max} / D_{\min}$ is lower than a critical value $ε *$ $ε^{*}$ , the solution of the second-order model is well approximated by the solution of the first-order Kuramoto model up to an error of the order of $ε$ $ε$ . They further argued that, for actual power systems, $ε$ $ε$ is of the order of 0.1 if one takes into account the control system in the effective damping constant $D i$ $D_{i}$ .

We now sketch the main idea of the proof of these theorems. For each point in time $t$ $t$ we identify the nodes $m$ $m$ and $n$ $n$ with the largest and smallest angles. That is, all angles $δ i (t)$ $δ_{i} (t)$ are in the arc $[δ n (t), δ m (t)]$ $[δ_{n} (t), δ_{m} (t)]$ , and

V (t) = δ m (t) - δ n (t)

$V (t) = δ_{m} (t) - δ_{n} (t)$

(129)

denotes the length of the arc. One can then compute the derivative $d V / d t$ $d V / d t$ explicitly and finds that this expression is strictly nonpositive if

Γ min sin (ζ) \geq cos (γ max) Γ crit .

$Γ_{\min} \sin (ζ) \geq \cos (γ_{\max}) Γ_{crit} .$

(130)

If this condition is satisfied, the arc $A (ζ)$ $A (ζ)$ cannot grow and the angles remain cohesive. Frequency synchronization is shown in a second step by differentiating Eq. (126) with respect to time as follows:

d d t ˙ δ i = - N \sum j = 1 K i j D i cos (δ i - δ j - γ i j) (˙ δ i - ˙ δ j) = - N \sum j = 1 ˜ Λ i j (t) δ j .

$\frac{d}{d t} {\dot{δ}}_{i} = - \sum_{j = 1}^{N} \frac{K_{i j}}{D_{i}} \cos (δ_{i} - δ_{j} - γ_{i j}) ({\dot{δ}}_{i} - {\dot{δ}}_{j}) = - \sum_{j = 1}^{N} {\tilde{Λ}}_{i j} (t) δ_{j} .$

(131)

Under the condition of phase cohesiveness, the matrix $˜ Λ (t)$ $\tilde{Λ} (t)$ is a time-varying directed Laplacian. Hence, Eq. (131) describes a contraction and all frequencies $˙ δ n$ ${\dot{δ}}_{n}$ must converge exponentially to a common value $Ω$ $Ω$ . Finally, the generalization to the second-order power system model is analyzed within singular perturbation theory.

2. Sufficient criteria for sparse networks

The previously presented sufficient criterion holds for all networks but is of little use if the topology is sparse since we then have $Γ min = 0$ $Γ_{\min} = 0$ . A different approach is needed in this case. One family of sufficient criteria for the existence of a phase-locked state uses a method introduced by Jadbabaie, Motee, and Barahona (2004). The steady-state condition (117) is reformulated in the following vectorial form:

P = ˜ Λ (δ) δ,

$P = \tilde{Λ} (δ) δ,$

(132)

where $˜ Λ (δ)$ $\tilde{Λ} (δ)$ is a state-dependent Laplacian with edge weights $w i j = K i j sin (δ i - δ j) / (δ i - δ j)$ $w_{i j} = K_{i j} \sin (δ_{i} - δ_{j}) / (δ_{i} - δ_{j})$ . One can then apply fixed point theorems to derive conditions guaranteeing that Eq. (132) has a solution. For instance, one can formulate the following condition in terms of the algebraic connectivity of the network (cf. Sec. II), which yields a sufficient condition for the existence of a synchronized state (Dörfler and Bullo, 2012a):

λ 2 (Λ) > ⎛ ⎝ \sum (i, j) \in E | P i - P j | 2 ⎞ ⎠ 1 / 2,

$λ_{2} (Λ) > {(\sum_{(i, j) \in E} | P_{i} - P_{j} |^{2})}^{1 / 2},$

(133)

where $λ 2 (Λ)$ $λ_{2} (Λ)$ is the algebraic connectivity of the ordinary weighted graph Laplacian. We note that one has to be careful about the domain when applying fixed point theorems. In particular, the sufficient conditions do not guarantee uniqueness of a solution, as discussed by Manik, Timme, and Witthaut (2017).

Another approach uses the cycle flow arguments that were introduced in Sec. III.A.6. Denoting the sine of the phase difference along the line $ℓ$ $ℓ$ as $ψ ℓ$ $ψ_{ℓ}$ , we have the general relation (36). It has been shown that the cycle flow contribution $f$ $f$ is actually negligible in most cases such that $ψ \approx E ⊤ Λ * P$ $ψ \approx E^{⊤} Λ^{*} P$ (Dörfler, Chertkov, and Bullo, 2013). Hence, a synchronous state with a phase cohesiveness $ζ$ $ζ$ exists if

sin (ζ) \geq ∥ ψ ∥ \infty \approx ∥ E ⊤ Λ * P ∥ \infty .

$\sin (ζ) \geq ∥ ψ ∥_{\infty} \approx ∥ E^{⊤} Λ^{*} P ∥_{\infty} .$

(134)

Equation (134) provides a rigorous sufficient condition for many elementary networks and an excellent approximate condition for actual power grid topologies (Dörfler, Chertkov, and Bullo, 2013). Recently generalized types of Lyapunov function were derived for lossless power grids (111), allowing researchers to derive sufficient conditions for the global stability of synchronous states (Schiffer, Efimov, and Ortega, 2018, 2019).

E. Probabilistic analysis of nonlinear stability

The synchronous state of a power grid can be globally stable in favorable cases, but typically several different asymptotic attractors exist. This includes limit cycles, as shown in Fig. 9 for a single generator. In this section we explore the nonlinear stability of different attractors of the swing equation using the probabilistic methods introduced in Sec. V.C.3. We note that the swing equation is a simplified model of power system dynamics in the vicinity of the synchronized state. Hence, nonsynchronous attractors of this model do not faithfully describe asymptotic states of real systems. However, the existence of these states shows that those aspects of the system well described by the swing equation are not sufficient to guarantee stability. This includes in particular inertia and primary droop control.

1. Limit cycles and other attractors

A large variety of only partially understood attractors exist in power system models. Limit cycles similar to the ones discussed for a single machine in Sec. V.C.2 are among the most common nonsynchronous attractors. These states can be understood as arising from a partial decoupling limit. Consider a network consisting of two disconnected parts $V l$ $V_{l}$ and $V r$ $V_{r}$ . The two parts can synchronize independently at respective frequencies $Ω l$ $Ω_{l}$ and $Ω r$ $Ω_{r}$ according to Eq. (116).

If the two parts are weakly coupled, the equations of motion include terms proportional to $sin [(Ω l - Ω r) t]$ $\sin [(Ω_{l} - Ω_{r}) t]$ , which average out in time. The limit cycle with separate frequencies $Ω l$ $Ω_{l}$ and $Ω l$ $Ω_{l}$ persists, but the coupling causes perturbations leading to oscillations around the mean frequencies (Gelbrecht, Kurths, and Hellmann, 2020). When the coupling increases, the basin size of these attractors (the variance of the frequency time series) may even increase until the attractors eventually become unstable. Figure 14 illustrates the nonlinear stability of such states in a network with two populations of oscillators with power injections $P i = \pm 1$ $P_{i} = \pm 1$ , effective damping $D i = 0.1$ $D_{i} = 0.1$ , and homogeneous line parameters $K$ $K$ . For weak coupling the oscillators always rotate freely at $˙ δ = P i / D i$ $\dot{δ} = P_{i} / D_{i}$ . When the coupling increases, a large variety of attractors come into being, most of them with negligible basin size (referred to as outliers). Two partially synchronous states of the type discussed earlier occur more often. The basin size first increases with $K$ $K$ , but the fully synchronized state then becomes dominant. Its basin size increases rapidly for $K \geq 4$ $K \geq 4$ , reaching unity for $K ≳ 8$ $K ≳ 8$ .

Nonlinear stability of different attractors in a small sample network. (a) Structure of the network with two classes of nodes with injections $P i = \pm 1$ $P_{i} = \pm 1$ . (b) The basin bifurcation diagram for the system with the synchronized, partially synchronized, and synchronized regimes. From Gelbrecht, Kurths, and Hellmann, 2020.

The occurrence of similar states was already studied from the point of view of hysteresis by Olmi et al. (2014). However, not all asymptotic states with an appreciable basin have this form. Olmi (2015) found that complex chimeras are possible, and Nitzbon et al. (2017) saw that perturbations at some types of nodes lead to limit cycles with a frequency not given by Eq. (116). The impact of noise, inertia, and network topology on different limit cycles and the hysteresis behavior was studied by Tumash, Olmi, and Schöll (2018).

2. The impact of network structure

We now return to the synchronized state and investigate its nonlinear stability with respect to localized perturbations in terms of the single-node basin stability. In particular, the set of initial states $A$ $A$ is chosen such that all nodes except one are at the fixed point values initially. This measure is of high practical relevance, as actual failures or disturbances are typically limited to a single machine or a small subnetwork. Further, it reveals the vulnerable nodes in a grid. A numerical analysis for sparse networks (Menck et al., 2014) showed that the local network features have a substantial impact on the single-node basin stability. In particular, dead ends tend to cause exceptionally poor single-node basin stability. Nodes at which the dead end is connected to the rest of the grid must be viewed as the most vulnerable spots with respect to dynamical perturbations.

The more detailed study of Schultz, Heitzig, and Kurths (2014a), based on a more realistic network ensemble introduced by Schultz, Heitzig, and Kurths (2014b), identified detours as stabilizing features and showed that there is a statistically robust relationship between network features and node robustness. Kim, Lee, and Holme (2015) showed that the community structure of the system determines which nodes reach a high single-node basin stability first when the coupling strength is increased. Kim, Lee, and Holme (2016) conducted a systematic exploration of the impact of network motifs on single-node basin stability and found that betweenness and degree matter specifically. Kim, Lee, and Holme (2016) and Kim et al. (2018) showed that there is strong nonmonotonic behavior of basin stability with increasing coupling strength. Finally, Kim et al. (2019) introduced integrated basin stability as a way to understand the transition to global synchrony as a function of increasing coupling strength.

A more detailed picture emerges if we add survivability, the probability that a perturbation does not lead to a violation of transient operational bounds. Hellmann et al. (2016) showed a strong direct dependence of single-node survivability on the node degree, with high degree nodes being particularly unstable. Combining survivability and basin stability, Nitzbon et al. (2017) showed a detailed picture of various single-node stabilities as a function of topology. Key to this is a node classification into those in a loop and those on a tree. Further subdividing the tree nodes into leaves, inner tree nodes, and sprouts of high and low neighbor degree led to the identification of classes of nodes that show qualitatively different stability behavior; see Fig. 15. Perturbations at some of these nodes were already shown to result in novel asymptotic states, mentioned in Sec. VI.E.1, that cannot be explained with a decoupling limit.

Topological properties determine the single-node basin stability and survivability. Left panel: survivability (Surv) and basin stability (BS) for nodes in a synthetic power grid network ensemble, showing pronounced differences in the behaviors of various node types. Right panel: node classification according to Nitzbon *et al.* (2017). Edges that are not in any cycles are trees. The nodes in trees are distinguished as roots (part of a cycle and a tree), sprouts (degree 1 nodes next to a root), leaf (degree 1 nodes not next to roots), and inner tree nodes. Sprouts are further distinguished according to their neighbor degree as dense ( $k a v > 5$ $k_{a v} > 5$ ) and sparse ( $k a v \leq 5$ $k_{a v} \leq 5$ ). Adapted from Nitzbon *et al.*, 2017.

3. The necessity of realistic models

When identifying the imprint of network structure on the stability properties of individual nodes, one usually proceeds by simplifying the system. To isolate the impact of topology, other factors are taken out. For example, most of the previously cited studies assumed homogeneous node parameters and neglected voltage dynamics and losses. We can use probabilistic methods to understand how much these assumptions actually alter the picture.

Wolff, Lind, and Maass (2018) showed that both inhomogeneity of the nodes and the details of the coupling of generators into the grid actually have a profound impact on the single-node stability. They evaluated the basin stability and the return times to the sync state for power grids with a fixed topology but while using a variety of models including the heterogeneity in the node and line parameters and the more detailed coupling of generators, where the generator was located at an internal node coupled by a line to the system bus. The grid was completely stable for homogeneous parameters but exhibited instability when introducing heterogeneity. The coupling of generators at internal nodes was found to stabilize the system overall; however, the location of unstable nodes also changed.

Taking higher-order internal dynamics of the generators into account for single-node stability analysis was considered by Auer et al. (2016). It was found that while single-node survivability is well captured by the swing equation (111), voltage dynamics and internal nodes can lead to additional asymptotic instabilities once outside the survivability region. The first results on basin stability of higher-order generator models were obtained by Liu et al. (2019).

Finally, most of the previously discussed works and much of the theoretical literature focus on lossless grids. Real transmission lines have a small but nonzero resistance (cf. Table I), and typical values of the parameter $γ$ $γ$ in Eq. (109) are around 0.24. While neglecting Ohmic losses is often justified for the power flow in the synchronized state, limit cycles and global stability properties are strongly affected (Hellmann et al., 2020). Figures 16(a)–16(d) show the basin of attraction of the synchronous state (peach) in the slice of phase space belonging to a single node centered on the sync state. As $γ$ $γ$ is increased from zero to realistic values, the basin first expands and then switches to the other half plane.

Main results of Hellmann *et al.* (2020). (a)–(d) Slices of the phase space in the dimensions belonging to a node $k$ $k$ . Different colors show different asymptotic states. Peach is the sync state, blue is a solitary state, and green is a solitary rotating against its power infeed. (e) Average single-node basin stability for the entire ensemble of networks, for the different types of asymptotic states, as well as the number of desynchronized oscillators following a perturbation (solid purple line). Adapted from Hellmann *et al.*, 2020.

A key role in this change is played by 1-solitary states, in which the entire network stays close to synchrony except for one desynchronized node. If the desynchronized node is rotating with the driving force $⟨ ˙ δ i ⟩ / P i > 0$ $⟨ {\dot{δ}}_{i} ⟩ / P_{i} > 0$ , we speak of a normal solitary. If it is rotating against the local driving force, we denote this as an exotic solitary state. Figure 16(e) shows the average asymptotic single-node basin stability for different types of asymptotic states: (i) peach for the sync state, (ii) dark blue for ordinary 1-solitary states, (iii) dark green for exotic 1-solitary states, (iv) light green and blue for more complex asymptotic states that also contain individual desynchronized nodes, and (v) gray for other asymptotic states.

The presence of exotic solitary states can be understood by considering the decoupling limit for an individual oscillator from the rest of the grid. When completely decoupled it will rotate with a frequency $Ω sol = P sol / D sol$ $Ω_{sol} = P_{sol} / D_{sol}$ . If we now increase the coupling, that state persists but, taking losses into account, the power flow no longer averages to 0 but instead to $P e = K e sin (γ e)$ $P_{e} = K_{e} \sin (γ_{e})$ according to Eq. (109). Hence, in contrast to the lossless case the coupling to the network will not only perturb the limit cycle but also shift its location to $Ω sol = [P sol - K e sin (γ)] / D sol$ $Ω_{sol} = [P_{sol} - K_{e} \sin (γ)] / D_{sol}$ . This picture was confirmed in detailed numerical experiments by Hellmann et al. (2020). Methods to restore synchronization from solitary states were discussed by Taher, Olmi, and Schöll (2019).

We stress that the aforementioned results were obtained for a particular power system model, and all models face certain limitations. For instance, the 1-solitary states are related to “pole slipping” transients in practice. They correspond to true limit cycles in the restricted model, but the dynamics is changed by the generator’s protection systems, which are not included in the model. We conclude that care should be taken when interpreting model simulation results for real-world situations. Models are useful in engineering power system stability. By construction, such models idealize certain real-world aspects and thus do not capture all details of real-world power system dynamics. In particular, many dynamical models do not account for the protection systems or changes in external conditions or control systems. One advantage of probabilistic methods for stability assessment is that they can be straightforwardly adapted to models of different complexity and scope.

VII. Structural Stability of Power Grids

Structural damages are the ultimate threat to power grid stability: “Typically, the blackout can be traced back to the outage of a single transmission (or generation) element,” according to Pourbeik, Kundur, and Taylor (2006). Such an initial outage can trigger secondary outages and, eventually, a cascade of failures bringing down a power grid entirely. One of the most important security regulations is the $N - 1$ $N - 1$ rule, with Wood, Wollenberg, and Sheblé, 2014 having stated that “no single outage will result in other components experiencing flow or voltage limit violations.”. But this rule is violated occasionally, making a grid vulnerable and large-scale blackouts possible.

In this section, we analyze the impact of structural damages focusing on outages of transmission elements. We mostly employ a quasistatic picture, assuming that after an outage the grid will rapidly relax to a new steady state (if it exists). But in this new state other transmission elements may be overloaded, leading to emergency shutdowns or short-circuit failures, etc.

The key questions we address in the following are (1) How do failures spread in the network? Given a transmission line failure, how and where are flows rerouted? (2) How does the entire network react to structural damages or other perturbations? Is there enough redundancy to cope with damages or do secondary failures take place? (3) Finally, how do large-scale blackouts emerge? How do cascades propagate through the grid?

A. Quasistatic analysis of line outages

1. The line outage distribution factors

We first analyze the outage of a single line within the linearized load flow approximation, or dc approximation, introduced in Sec. III.A.4. Assume that a line $ℓ^= (r, s)$ $ℓ \hat{=} (r, s)$ that initially carries the flow $F (0) ℓ$ $F_{ℓ}^{(0)}$ fails. If the grid remains connected, the flow change at another transmission line $e^= (m, n)$ $e \hat{=} (m, n)$ is linear in $F (0) ℓ$ $F_{ℓ}^{(0)}$ ,

Δ F e = LODF e, ℓ F (0) ℓ,

$Δ F_{e} = {LODF}_{e, ℓ} F_{ℓ}^{(0)},$

where the factor of proportionality is referred to as a line outage distribution factor (LODF) (Grainger and Stevenson, 1994; Wood, Wollenberg, and Sheblé, 2014).

We summarize the derivation of the LODFs following Guo et al. (2009) and Strake et al. (2019). As the line $ℓ = (r, s)$ $ℓ = (r, s)$ fails, the nodal susceptance matrix changes as

B \to B + Δ B = B + b r s ν r s ν r s ⊤,

$B \to B + Δ B = B + b_{r s} ν_{r s} {ν_{r s}}^{⊤},$

(135)

where the vector $ν r s \in R N$ $ν_{r s} \in R^{N}$ is $+ 1$ $+ 1$ at position $r$ $r$ , $- 1$ $- 1$ at position $s$ $s$ , and 0 otherwise. This causes a change of the nodal phase $θ \to θ' = θ + α$ $θ \to θ^{'} = θ + α$ . Subtracting the linearized load flow equations (32) for the perturbed and unperturbed grid then yields

(B + Δ B) α = - Δ B θ,

$(B + Δ B) α = - Δ B θ,$

(136)

which can be solved for $α$ $α$ and used to compute the change of the line flows as

Δ F = B d E ⊤ α = - B d E ⊤ (B + Δ B) * Δ B θ,

$Δ F = B_{d} E^{⊤} α = - B_{d} E^{⊤} (B + Δ B)^{*} Δ B θ,$

(137)

with the superscript $*$ $*$ denoting the Moore-Penrose pseudoinverse. Equation (137) can be greatly simplified using the Woodbury matrix identity, which finally yields

Δ F = B d E ⊤ B * ν r s (1 - b r s ν r s ⊤ B * ν r s) - 1 F (0) r s,

$Δ F = B_{d} E^{⊤} B^{*} ν_{r s} (1 - b_{r s} {ν_{r s}}^{⊤} B^{*} ν_{r s})^{- 1} F_{r s}^{(0)},$

(138)

and thus $LODF e, ℓ = b m n ν ⊤ m n B * ν r s (1 - b r s ν r s ⊤ B * ν r s) - 1$ ${LODF}_{e, ℓ} = b_{m n} ν_{m n}^{⊤} B^{*} ν_{r s} (1 - b_{r s} {ν_{r s}}^{⊤} B^{*} ν_{r s})^{- 1}$ .

This procedure is readily generalized to multiple line outages (Güler, Gross, and Liu, 2007; Kaiser, Strake, and Witthaut, 2020). Assume that the lines $ℓ 1, ℓ 2, \dots, ℓ M$ $ℓ_{1}, ℓ_{2}, \dots, ℓ_{M}$ fail, but the network remains connected. We define a projection matrix from the space of all links onto the subset of failing links $P \in R M \times L$ $P \in R^{M \times L}$ as follows:

P m e = {1 if e = ℓ m, 0 otherwise .

$P_{m e} = {\begin{cases} 1 & if e = ℓ_{m}, \\ 0 & otherwise . \end{cases}$

We define projections of the node-edge incidence matrix, the branch reactance matrix, and initial flow vectors as

D ≔ P E \in R N \times M, B out ≔ P B d P ⊤ = diag (b ℓ 1, b ℓ 2, \dots, b ℓ M) \in R M \times M, F (0) out = P F (0) = (F (0) ℓ 1, \dots, F (0) ℓ M) \in R M .

$D ≔ P E \in R^{N \times M}, B_{out} ≔ P B_{d} P^{⊤} = diag (b_{ℓ_{1}}, b_{ℓ_{2}}, \dots, b_{ℓ_{M}}) \in R^{M \times M}, F_{out}^{(0)} = P F^{(0)} = (F_{ℓ_{1}}^{(0)}, \dots, F_{ℓ_{M}}^{(0)}) \in R^{M} .$

One can then proceed as for a single failing line. Using the Woodbury matrix identity, one obtains

Δ F = B d E ⊤ B * D (1 - B out D ⊤ B * D) * F (0) out .

$Δ F = B_{d} E^{⊤} B^{*} D {(1 - B_{out} D^{⊤} B^{*} D)}^{*} F_{out}^{(0)} .$

(139)

In Eq. (139), we only have to invert an $M \times M$ $M \times M$ matrix in addition to the inversion of the initial matrix $B *$ $B^{*}$ . If $M$ $M$ is small, the inverse can be evaluated explicitly and we obtain a set of generalized LODFs.

We note that contingency analysis via LODFs can be improved by a modification of the linearization procedure. One starts with the nonlinear expression $F j \to k = b j k sin (θ j - θ k)$ $F_{j \to k} = b_{j k} \sin (θ_{j} - θ_{k})$ for the real power flow and carries out the linearization only at a later stage (Jung and Kettemann, 2016; Manik et al., 2017). One then obtains the same results as given previously, but with reweighted line susceptances

b j k \to b j k cos (θ j - θ k) .

$b_{j k} \to b_{j k} \cos (θ_{j} - θ_{k}) .$

(140)

Related approaches are sometimes referred to as a “hot-start dc approximation.”

2. Spreading of failures

A deeper physical insight into the network flow rerouting problem is obtained by an analogy to discrete electrostatics. Using again the Woodbury matrix identity, Eq. (136) can be recast in the form

B α = q,

$B α = q,$

(141)

with the source term

q = (1 - b r s ν r s ⊤ B * ν r s      ≕ χ r s) - 1 F (0) r s ν s r .

$q = (1 - \underset{≕ χ_{r s}}{\underset{⏟}{b_{r s} {ν_{r s}}^{⊤} B^{*} ν_{r s}}})^{- 1} F_{r s}^{(0)} ν_{s r} .$

(142)

As the matrix $B$ $B$ is a Laplacian matrix and the right-hand side is nonzero only at positions $r$ $r$ and $s$ $s$ with opposite sign, Eq. (141) is a discrete Poisson equation with a dipole source and $α$ $α$ is a dipole potential; see Biggs (1997). Indeed, if we compute the impact of the outage of a line in a rectangular square grid, we simply recover a discretized version of the familiar electric dipole field, as illustrated in Fig. 17. Thus, $LODF e, ℓ$ ${LODF}_{e, ℓ}$ depends on the orientation of the two lines $e$ $e$ and $ℓ$ $ℓ$ and decreases as $(distance) - 2$ $({distance)}^{- 2}$ .

Impact of a link failure in a square lattice with uniform edge weights. (a) Normalized change of the nodal potentials $α n$ $α_{n}$ for a single failing link located in the center of the network. Both the size of the nodes and the color code represent $α n$ $α_{n}$ . (b) Normalized change of the link flows $Δ F e$ ${Δ F}_{e}$ for the same topology. Arrows and color represent the direction and strength of flow changes, respectively. From Strake *et al.*, 2019.

Understanding the impact of line outages for arbitrary grid topologies is much more challenging. The main complexity arises from the network topology encoded in the Laplacian $B$ $B$ , which can be highly irregular. Still we can get some insight using tools from algebraic graph theory (Strake et al., 2019; Kaiser and Witthaut, 2021a).

Before we proceed to spatial failure spreading, we have a more detailed look at the right-hand side of the Poisson equation (141). The strength of the dipole $| q |$ $| q |$ , or more precisely the factor $χ r s$ $χ_{r s}$ , can be understood from the grid topology: If a unit power is injected at node $r$ $r$ and withdrawn at node $s$ $s$ , then $χ r s$ $χ_{r s}$ is the flow via the direct link $(r, s)$ $(r, s)$ and $1 - χ r s$ $1 - χ_{r s}$ is the flow via other nondirect pathways. Hence, we can view $1 - χ r s$ $1 - χ_{r s}$ as a measure of redundancy of the line $(r, s)$ $(r, s)$ . Indeed, it can be rigorously related to topological redundancy measures (Strake et al., 2019; Guo et al., 2021). An alternative approach to line outage problems was introduced by Ronellenfitsch et al. (2017). Flow changes due to a failure can be decomposed into cycle flows as $Δ F = C f$ $Δ F = C f$ using the edge-cycle incidence matrix (5). The cycle flow strength $f$ $f$ is again determined by a discrete Poisson equation, now formulated on the dual graph.

3. Locality and the importance of distance

Intuitively, a transmission line outage should affect nearby lines more heavily than remote ones. In a regular square lattice, flow changes decay as $(distance) - 2$ ${(distance)}^{- 2}$ , as previously discussed. But how can we quantify this statement and what does nearby mean in a network with complex topology? These questions were studied using a spectral approach by Labavic et al. (2014), Jung and Kettemann (2016), Kettemann (2016), and Rohden et al. (2016).

The discrete Poisson equation (141) can be formally solved by decomposing the Laplacian $B$ $B$ into its eigenvectors defined as $B Φ j = λ j Φ j$ $B Φ_{j} = λ_{j} Φ_{j}$ with indices $j = 1, \dots, N$ $j = 1, \dots, N$ . Assuming the failure of a line $(r, s)$ $(r, s)$ , we obtain

α m = - Δ b r s sin (θ r - θ s) N \sum j = 2 Φ j m (Φ * j r - Φ * j s) λ j

$α_{m} = - Δ b_{r s} \sin (θ_{r} - θ_{s}) \sum_{j = 2}^{N} \frac{Φ_{j m} (Φ_{j r}^{*} - Φ_{j s}^{*})}{λ_{j}}$

(143)

using the generalized hot-start linear response with edge weights (140); see also Haehne et al. (2019). Similarly, the flow changes at the remaining links are given by

Δ F m n = b m n cos (θ m - θ n) Δ b r s sin (θ r - θ s) \times N \sum j = 2 (Φ j m - Φ j n) (Φ * j r - Φ * j s) λ j .

$Δ F_{m n} = b_{m n} \cos (θ_{m} - θ_{n}) Δ b_{r s} \sin (θ_{r} - θ_{s}) \times \sum_{j = 2}^{N} \frac{(Φ_{j m} - Φ_{j n}) (Φ_{j r}^{*} - Φ_{j s}^{*})}{λ_{j}} .$

(144)

For a regular grid of degree $2 d$ $2 d$ , the eigenvectors $Φ j$ $Φ_{j}$ are plain waves and the eigenvalues $λ j$ $λ_{j}$ can be obtained exactly, yielding (Kettemann, 2016)

Δ F m n \sim A m n Δ b r s c d | m - r | - d .

$Δ F_{m n} \sim A_{m n} Δ b_{r s} c_{d} | m - r |^{- d} .$

(145)

In Eq. (145) $c d$ $c_{d}$ are constants. For $d = 2$ $d = 2$ , a regular square grid, the power law decays with exponent 2 and $c 2 = 2 π$ $c_{2} = 2 π$ . For a random graph one finds exponentially localized eigenstates, and thus an exponentially decaying response (Torres-Sánchez, Freitas de Abreu, and Kettemann, 2020), so the disturbance due to the outage affects the power grid in only a finite range around the outage (Kettemann, 2016). This phenomenon is well known and originates from random scattering from nodes with random degree (Hata and Nakao, 2017) and/or random parameters (García-Mata et al., 2017), so-called Anderson localization (Anderson, 1958).

Numerical calculations based on the ac load flow equations were found to be in good agreement with Eq. (143); see Jung and Kettemann (2016). That study considered the addition of a new transmission line and confirmed the generic power-law decay of the change in transmitted power $Δ F \sim (distance) - 2$ $Δ F \sim {(distance)}^{- 2}$ up to finite size corrections. As a more realistic grid topology, the effect of a change in transmission line capacity was studied in a model of the German transmission grid (Medjroubi, Matke, and Kleinhans, 2015) that is shown in Fig. 18 (left image). The largest 2-connected component of this grid was considered, with $N = 260$ $N = 260$ nodes and $L = 479$ $L = 479$ edges, homogeneous line parameter $b i j$ $b_{i j}$ and power injections randomly chosen from a binary distribution, $P i \in {- P, P}$ $P_{i} \in {- P, P}$ . Distance is defined as the geodesic distance; see Sec. II. The resulting response behavior depends strongly on the location of the added line. Most locations result in a long-range power-law decay with an exponent that is, on average, close to 2. However, there are certain classes of locations, denoted in Fig. 18 as subsets 1 and 2, where the change results in an exponential decay beyond a certain length scale. Notably, the regions around subset 2 are weakly connected to the remaining grid, forming weakly connected islands, which explains the fact that the disturbance decays exponentially at distances beyond these islands.

Left panel: German transmission grid topology (220 and 380 kV) (Medjroubi, Matke, and Kleinhans, 2015). Right panel: Double-logarithmic plot of transmitted power change $⟨ | Δ F (ℓ) m n | ⟩ (r)$ $⟨ | Δ F_{m n}^{(ℓ)} | ⟩ (r)$ as a function of distance $r$ $r$ when a transmission line $(ℓ)$ $(ℓ)$ between 880 different pairs of nodes $i$ $i$ and $j$ $j$ was added. We averaged this quantity over all edges with the same distance $r$ $r$ to the added line $(ℓ)$ $(ℓ)$ and performed an ensemble average over $R = 100$ $R = 100$ realizations with randomly distributed generators and loads. The power-law expected for a two-dimensional regular grid $Δ F \sim r - 2$ $Δ F \sim r^{- 2}$ is shown for comparison (black line). The colored data belong to special subsets of added edges, as discussed in the text. Error bars, 95% confidence level. Adapted from Jung and Kettemann, 2016.

A different approach to understanding the spatial pattern of flow rerouting is to adapt the measure of distance to the respective problem (Strake et al., 2019; Kaiser and Witthaut, 2021a). Consider the line outage distribution factors $LODF e, ℓ$ ${LODF}_{e, ℓ}$ for the IEEE 30-node test grid shown in Fig. 19, fixing a failing link $ℓ$ $ℓ$ in the right part of the network. The LODFs decay with the geodesic distance of two links, but the correlation is only moderately pronounced. This is not surprising, as the geodesic distance does not take into account the fact that flow rerouting requires two paths between the links $ℓ$ $ℓ$ and $e$ $e$ , one to and one fro. We can define an alternative distance measure, the rerouting distance, that seeks the shortest path from one end of the link $ℓ$ $ℓ$ via $e$ $e$ to the other end of the link $ℓ$ $ℓ$ . Numerical studies show that the correlation of rerouting distance to the LODFs is much stronger than for the geodesic distance (Strake et al., 2019).

Localization of the impact of line failures as seen with different distance measures. (a) Flow changes after the outage of a single line (in red) in the IEEE 30-bus test grid. (b) Geodesic distance to the failing link. (c) Rerouting distance to the failing link accurately predicts the spatial distribution of the flow changes. (d),(e) Flow changes $| Δ F e |$ $| Δ F_{e} |$ vs distance to the failing link $ℓ$ $ℓ$ on a logarithmic scale. One clearly observes that the rerouting distance used in (e) yields a higher correlation than the ordinary geodesic distance shown in (d). Adapted from Strake *et al.*, 2019.

Finding a perfect measure of distance requires one to take all details of the Laplacian $B$ $B$ into account, not just the single shortest path. Indeed, the resistance distance (Klein and Randić, 1993) provides such a measure of distance and can be applied to line outages (Tyloo, Pagnier, and Jacquod, 2019). However, as this distance is calculated from the full matrix $B$ $B$ , it does not offer a direct topological interpretation as the shortest path distance does.

4. The importance of community structures

Some power grids show a pronounced community structure. In the Scandinavian grid, for example, Finland is connected to its neighbors via only two ac lines [Fig. 20(a)]. This weak topological connectivity also induces a weak electric connectivity: LODFs between two lines in different communities are significantly smaller than expected from their distance. This aspect is exemplarily illustrated for the failure of a line in Sweden [Fig. 20(b)]. For most lines in Finland the respective LODFs are of the order of $10 - 7$ $10^{- 7}$ , but some lines in the north near the connection to Sweden are slightly more affected. A more quantitative analysis is provided in Fig. 20(c) showing that LODFs between different communities are significantly suppressed, up to several orders of magnitude, compared to LODFs within a community. Hence, a splitting into communities can strongly impede the spreading of failures (Kaiser et al., 2021). Several other topological structures have similar or even stronger effects on failure spreading (Kaiser, Latora, and Witthaut, 2021; Kaiser and Witthaut, 2021a): If the coupling between two parts of a grid can be described using an adjacency matrix of unit rank, then all mutual LODFs vanish exactly and the spreading of failures and cascades is suppressed entirely. The simplest realization of such a unit-rank coupling is given by a bridge (Ronellenfitsch et al., 2017; Guo et al., 2021) that is exploited in recent proposals for emergency measures to contain cascading failures (Zocca et al., 2021; Bialek and Vahidinasab, 2022). These results emphasize the importance of structural features for the general robustness of electric power grids.

The importance of community structures in the spreading of failures. (a) Community structure of the Scandinavian power grid obtained by spectral clustering. (b) Magnitude of the LODFs for a failing link at the border between southern Norway and Sweden (red line). (c) Magnitude of the line outage distribution factors $| LODF j ℓ |$ $| {LODF}_{j ℓ} |$ as a function of the unweighted edge distance of the transmission lines $j, ℓ$ $j, ℓ$ . The cases in which $j$ $j$ and $ℓ$ $ℓ$ are in the same or different communities are compared. The line gives the median, while the shaded area gives the 25%–75% quantile. The grid data were taken from Hörsch, Hofmann *et al.* (2018).

5. Multiple outages and collective effects

When two lines $k$ $k$ and $ℓ$ $ℓ$ fail, we cannot simply superpose the flow changes but instead have to account for collective effects. To see this, assume that $k$ $k$ and $ℓ$ $ℓ$ fail successively. In calculating the impact of the failure $ℓ$ $ℓ$ , we must take into account the observation that the flow $F ℓ$ $F_{ℓ}$ has already changed due to the failure of $k$ $k$ . The correct result for the impact of a multiple line outage was derived in Sec. VII.A.1 and, for two lines, reduces to

Δ F (k, ℓ) e = (L e ℓ L e k) (1 - L ℓ k - L k ℓ 1) - 1 ⎛ ⎝ F (0) ℓ F (0) k ⎞ ⎠,

$Δ F_{e}^{(k, ℓ)} = (\begin{matrix} L_{e ℓ} & L_{e k} \end{matrix}) {(\begin{matrix} 1 & - L_{ℓ k} \\ - L_{k ℓ} & 1 \end{matrix})}^{- 1} (\begin{matrix} F_{ℓ}^{(0)} \\ F_{k}^{(0)} \end{matrix}),$

(146)

using the shorthand $L e ℓ = LODF e, ℓ$ $L_{e ℓ} = {LODF}_{e, ℓ}$ and denoting the failing lines in the superscript. Collective effects can have surprising consequences. The first impact of a line failure may be partly compensated for by the intentional removal of a second line. Thus, suppose that the failure of line $ℓ$ $ℓ$ has led to an increase of the flow on a critical line $e$ $e$ . In many cases we can find another line $k$ $k$ such that

| F (k, ℓ) e | \leq | F (ℓ) e | .

$| F_{e}^{(k, ℓ)} | \leq | F_{e}^{(ℓ)} | .$

(147)

This concept was previously discussed by Motter (2004) and Witthaut and Timme (2015) for different types of flow networks; we return to this in Sec. VII.D. More surprisingly, we can find cases where the collective impact of two outages is the opposite of the impact of two isolated outages, i.e.,

Δ F (k) e, Δ F (ℓ) e > 0 but Δ F (k, ℓ) e < 0 .

$Δ F_{e}^{(k)}, Δ F_{e}^{(ℓ)} > 0 but Δ F_{e}^{(k, ℓ)} < 0 .$

(148)

An example of this effect is shown in Fig. 21, and a more detailed discussion was given by Kaiser, Strake, and Witthaut (2020). Further analysis shows that collective effects are particularly strong whenever the mutual LODFs, or more precisely the expression $\sqrt L k ℓ L ℓ k$ $\sqrt{L_{k ℓ} L_{ℓ k}}$ , are large.

Collective effects can lead to a complete reversal of the flow changes in an $N - 2$ $N - 2$ failure. The color code of the lines indicates the magnitude of flow, with red indicating failing links. (a) Flows in the initial unperturbed grid. (b),(c) Flow changes after the individual failure of two links ( $ℓ$ $ℓ$ and $k$ $k$ , respectively, in red). In both cases, the flow in the top left link ( $e$ $e$ ; boldface font) is greater than in the unperturbed grid: $Δ F (ℓ) l \approx 0.3$ $Δ F_{l}^{(ℓ)} \approx 0.3$ , $Δ F (k) l \approx 0.1$ $Δ F_{l}^{(k)} \approx 0.1$ . (d) Flows after the simultaneous failure of both links $ℓ$ $ℓ$ and $k$ $k$ . The flow in the top left link is smaller in magnitude than in the unperturbed grid and even changes its direction: $Δ F (ℓ, k) e \approx - 0.7$ $Δ F_{e}^{(ℓ, k)} \approx - 0.7$ . From Kaiser, Strake, and Witthaut, 2020.

B. Robustness of power grids and critical infrastructures

The outage of a transmission element leads to the rerouting of currents and power flows, as analyzed in Sec. VII.A. This may eventually cause secondary failures of other transmission elements, which can result in a large-scale cascade affecting the entire grid. In this section we investigate the grid vulnerability to secondary failures and analyze the weak spots of a network.

1. Why secondary failures?

If the current on an overhead transmission line is too large, it will heat up due to Ohmic losses, which may lead to line sag and eventually a short-circuit fault, as discussed in Sec. V.B. Strict security rules are commonly implemented to avoid such faults such that the current on a line $ℓ$ $ℓ$ may not exceed a certain limit, the line rating $| I ℓ | \leq I max ℓ$ $| I_{ℓ} | \leq I_{ℓ}^{\max}$ . If voltages are close to the set point and losses are small, this directly translates into an upper limit for the real power flow

| F ℓ | \leq F max ℓ .

$| F_{ℓ} | \leq F_{ℓ}^{\max} .$

(149)

A violation of this limit typically leads to an emergency shutdown if the line is appropriately monitored. The line will not be damaged, but it is no longer available for power transmission. In practice, one differentiates between short-term and long-term loads of a transmission line. A higher loading may be acceptable on short terms, such that a higher line rating applies. Furthermore, not only currents but also voltages are important for grid stability. Strict security rules apply to the voltage level at every substation, which must remain within a certain interval around the grid reference level. A decreasing voltage may indicate a looming voltage instability, as discussed in Sec. V.A.

2. Critical links: A graph theoretic perspective

Consider a heavily loaded power transmission grid which is no longer $N - 1$ $N - 1$ secure. Some transmission lines may be critical in the sense that if they fail secondary overloads occur. But which lines are critical and which lines are prone to secondary overloads? These questions can be answered with extensive numerical simulations, but a better analytic understanding would be helpful. We expect outages of lines with a high load or flow to be more harmful, but this is certainly not the entire story. Whether the grid bears enough redundancy to cope with the failure is equally important. These aspects can be quantitatively understood using tools from graph theory (Witthaut et al., 2016).

Assume that a line $(a, b)$ $(a, b)$ fails that initially carried the real-power flow $F (0) a \to b$ $F_{a \to b}^{(0)}$ from node $a$ $a$ to node $b$ $b$ . After the failure, the power must be rerouted from $a$ $a$ to $b$ $b$ via different pathways. Graph theory now provides a necessary condition for this to be possible. The Edmonds-Karp algorithm yields the maximum flow that can be transmitted from $a$ $a$ to $b$ $b$ while respecting the line limits, which we call the redundant capacity $K red a \to b$ $K_{a \to b}^{red}$ of the line $(a, b)$ $(a, b)$ . The algorithm also yields the set of edges limiting the flow from $a$ $a$ to $b$ $b$ , which is referred to as the minimum cut in graph theory. Hence, we can obtain both a measure of redundancy and the locations of the bottlenecks in the grid.

In summary, a necessary condition for having a feasible power flow after the outage of line $(a, b)$ $(a, b)$ is given by $| F (0) a \to b / K red a \to b | \leq 1$ $| F_{a \to b}^{(0)} / K_{a \to b}^{red} | \leq 1$ . In practice, secondary overloads will typically occur earlier, as power flows are generally not maximum flow network flows but have to respect Kirchhoff’s laws. An example is shown in Fig. 22 for a semisynthetic power grid based on the structure of the Austrian transmission grid. One of the marked links carries a high real power flow of about 1.2 GW. It is uncritical because it has a high redundant capacity. The second line has a lower flow ( $\leq 0.9 GW$ $\leq 0.9 GW$ ) but lacks redundancy. The ratio $| F (0) a \to b / K red a \to b |$ $| F_{a \to b}^{(0)} / K_{a \to b}^{red} |$ is larger than 1 such that an outage of this line will always lead to secondary outages and a fragmentation of the grid. Applying the Edmonds-Karp algorithm directly yields the bottleneck that limits the redundancy.

(a) Real power flows $F a b$ $F_{a b}$ in MW in a test grid based on the Austrian grid topology. The two marked links are investigated in detail. (b) The ratio of the flow and the redundant capacity $F (0) a \to b / K red a \to b$ $F_{a \to b}^{(0)} / K_{a \to b}^{red}$ provides a powerful indicator for critical links. (c),(d) Illustration of rerouting pathways after the outage of links 1 and 2, respectively. Several rerouting pathways exist for link 1 such that the redundant capacity $K red$ $K^{red}$ is rather high, and a line failure is uncritical. In contrast, a failure of link 2 will cause secondary outages. The min cut for the rerouting problem consists of a single link (marked with an arrow). This bottleneck has a free capacity of 864 MW such that a rerouting is impossible. Grid data are based on the work of Hutcheon and Bialek (2013), with a slight adaptation of power injections. Austria was artificially islanded for illustrative purposes.

Extensive numerical tests for different power system models (Witthaut et al., 2016) show that the ratio $| F (0) a \to b / K red a \to b |$ $| F_{a \to b}^{(0)} / K_{a \to b}^{red} |$ is a powerful predictor for critical links. Furthermore, one can improve power system robustness by strengthening the bottlenecks (Rohden et al., 2017).

3. Generation variability and critical fluctuations

Thus far we have investigated potential overloads and cascades assuming that the power injections are fixed and perfectly known. But even if we know the scheduled values, there are always some residual fluctuations of generation and demand. How do these fluctuations affect power system stability? Is it possible that the grid is safe on average but that overloads occur during stochastic peaks that may eventually trigger a cascade?

These important questions were addressed by Nesti, Zocca, and Zwart (2018) using a quasistatic picture. The power injections $P$ $P$ are treated as Gaussian random variables with a mean $μ P$ $μ_{P}$ and a correlation matrix $ε Σ P$ $ε Σ_{P}$ . The essential quantities in a security assessment are the relative line loads $L ℓ ≔ F ℓ / F max ℓ$ $L_{ℓ} ≔ F_{ℓ} / F_{ℓ}^{\max}$ , which are determined by the power injections via

L = diag (F max 1, \dots, F max L) - 1 B d E ⊤ B *      ≕ Υ P

$L = \underset{≕ ϒ}{\underset{⏟}{diag (F_{1}^{\max}, \dots, F_{L}^{\max})^{- 1} B_{d} E^{⊤} B^{*}}} P$

using Eq. (33). Hence, the line loads $L$ $L$ are also Gaussian random variables with mean $μ L = Υ μ P$ $μ_{L} = ϒ μ_{P}$ and correlation matrix $ε Σ L = ε Υ Σ P Υ ⊤$ $ε Σ_{L} = ε ϒ Σ_{P} ϒ^{⊤}$ . A line $ℓ$ $ℓ$ is overloaded if $L ℓ > 1$ $L_{ℓ} > 1$ . In the limit of small fluctuations $ε \to 0$ $ε \to 0$ one can derive the following probability for this event using large-deviation theory:

P (| L ℓ | \geq 1) \approx exp (- (1 - | μ L, ℓ |) 2 2 ε σ 2 ℓ),

$P (| L_{ℓ} | \geq 1) \approx \exp (- \frac{(1 - | μ_{L, ℓ} |)^{2}}{2 ε σ_{ℓ}^{2}}),$

(150)

where $σ 2 ℓ = (Σ L) ℓ ℓ$ $σ_{ℓ}^{2} = (Σ_{L})_{ℓ ℓ}$ . While this is an approximation for nonvanishing fluctuations, it was found useful for ranking purposes. The smaller the ratio $(1 - | μ L, ℓ |) 2 / σ 2 ℓ$ $(1 - | μ_{L, ℓ} |)^{2} / σ_{ℓ}^{2}$ , the more likely an overload of a line $ℓ$ $ℓ$ is.

An example of such a probabilistic contingency analysis is shown in Fig. 23 for the SciGRID model of the German power grid. The nominal power injections were determined via an OPF, and the correlation matrix $Σ P$ $Σ_{P}$ was derived from a statistic analysis of actual generation time series. One finds that the vulnerable lines where $P (| L ℓ | \geq 1)$ $P (| L_{ℓ} | \geq 1)$ is largest do not coincide with the most heavily loaded lines. High loads are one important indicator for vulnerability, but the network structure and the generation variability enter the overload probability on an equal footing via the effective variance $σ 2 ℓ$ $σ_{ℓ}^{2}$ .

Fluctuations of power generation and consumption can induce overloads of transmission lines. Left panel: nominal line loads $| μ L, ℓ | = | ⟨ L ℓ ⟩ |$ $| μ_{L, ℓ} | = | ⟨ L_{ℓ} ⟩ |$ in the SciGRID model of the German power grid. The nominal operation is determined from one-hour averages of the renewable feed-in and load. The feed-in of the dispatchable power plants and the power flows are then computed via a linear OPF; see Sec. III.B. Right panel: probability of a transmission line overload $P (| L ℓ | \geq 1)$ $P (| L_{ℓ} | \geq 1)$ when the true power injections fluctuate around the nominal values. Vulnerable lines are identified using small values of the indicator $(1 - | μ L, ℓ |) 2 / σ 2 ℓ$ $(1 - | μ_{L, ℓ} |)^{2} / σ_{ℓ}^{2}$ . From Nesti, Zocca, and Zwart, 2018.

One can extend this analysis and ask which power injections typically lead to the failure of a single line. In mathematical terms: What is the conditional probability distribution of $P$ $P$ given that $| L ℓ | \geq 1$ $| L_{ℓ} | \geq 1$ , i.e., given that the line $ℓ$ $ℓ$ is overloaded? For $ε \to 0$ $ε \to 0$ this probability distribution gets sharply concentrated around the vector

P (ℓ) = μ P + σ - 2 ℓ [sgn (μ L, ℓ) - μ L, ℓ] Σ P Υ ⊤ e ℓ,

$P^{(ℓ)} = μ_{P} + σ_{ℓ}^{- 2} [sgn (μ_{L, ℓ}) - μ_{L, ℓ}] Σ_{P} ϒ^{⊤} e_{ℓ},$

(151)

where sgn denotes the sign function and $e ℓ$ $e_{ℓ}$ is the $ℓ$ $ℓ$ th unit vector. That is, there is a typical injection pattern leading to an overload of a line given that fluctuations are small but nonzero (which is a plausible assumption for large power transmission grids). An important conclusion is that a line failure is not necessarily triggered by large fluctuations in the neighborhood but rather can be triggered by the cumulative effects of small unusual fluctuations in the entire network. Correlations of the power injections, such as by a sudden large-scale increase in wind power generation, play a key role in such an outage.

C. Cascades of failures and large-scale blackouts

The outage of a transmission line leads to a rerouting of power flows that may cause secondary overloads, as previously discussed. In this section we discuss how this mechanism leads to a cascade of failures and how such a cascade propagates through the grid.

1. Simulating cascading failures

Cascading failures are often simulated in a quasistatic picture. One computes a series of stationary states of the grid by solving the ac [Eqs. (25)] or dc [Eqs. (32)] load flow equations, where the topology changes as elements of the grid fail. The basic steps of a cascade may be simulated using Algorithm 1.

Algorithm 1.

for all branches

t

$t$ that may fail do

Remove trigger link

t

$t$

repeat

Solve ac or dc load flow equations

for all branches

ℓ

$ℓ$ do

if branch

ℓ

$ℓ$ is overloaded then

Remove branch

ℓ

$ℓ$ from grid

end if

end for

until No further overloads or grid disconnected

end for

More advanced simulations can include aspects of voltage stability as well as actions by the grid operators such as load shedding, the immediate disconnection of consumers to reduce the grid load.

We note that a wide body of literature on cascading failures exists in the statistical physics literature. Purely topological approaches based on percolation theory are particularly popular, as they allow for analytic solutions; see Albert, Jeong, and Barabási (2000), Albert, Albert, and Nakarado (2004), and Newman (2012). However, the applicability of topological models to power grid stability is limited, as discussed by Hines, Cotilla-Sanchez, and Blumsack (2010), Bompard, Luo, and Pons (2015), and Korkali et al. (2017). Topological models and measures can provide some general insight into the vulnerability of networks [see Galindo-González, Angulo-García, and Osorio (2020)], but they may be misleading when applied to particular grids. In particular, they capture neither the physics of power flows nor the heterogeneity and localization of power generation and consumption.

2. Local versus nonlocal propagation

How do cascading failures propagate through a power grid? Do subsequent failures occur in close proximity, or can they jump to remote areas of a grid?

A prime example of a locally propagating cascade was observed during the 2006 western European blackout (Union for the Coordination of Transmission of Electricity, 2007). The cascade was triggered by the shutdown of a line and a switching event in northwestern Germany. In every step of the cascade, flow was rerouted to the next available route in the southeastern direction, thereby causing secondary failures along these routes.

However, secondary outages can also take place at rather long distances from the initial failure. A simple reason is that small flow changes can also cause a line outage if a line were already heavily loaded before. In particular, a failure of line $ℓ$ $ℓ$ causes a secondary overload of line $e$ $e$ if (assuming that $F (0) e > 0$ $F_{e}^{(0)} > 0$ without loss of generality)

F (0) e + LODF e, ℓ F (0) ℓ > F max e \Leftrightarrow LODF e, ℓ F (0) ℓ > F max e - F (0) e .

$F_{e}^{(0)} + {LODF}_{e, ℓ} F_{ℓ}^{(0)} > F_{e}^{\max} \Leftrightarrow {LODF}_{e, ℓ} F_{ℓ}^{(0)} > F_{e}^{\max} - F_{e}^{(0)} .$

(152)

Hence, a secondary outage occurs if either the flow change $Δ F e = LODF e, ℓ F (0) ℓ$ $Δ F_{e} = {LODF}_{e, ℓ} F_{ℓ}^{(0)}$ is large or the line $e$ $e$ was heavily loaded before the outage such that $F max e - F (0) e$ $F_{e}^{\max} - F_{e}^{(0)}$ is small. While the flow changes $Δ F e$ $Δ F_{e}$ are strongest in the vicinity of the failing line $ℓ$ $ℓ$ , little can be said about the initial line loadings $F (0) e$ $F_{e}^{(0)}$ . They are determined by the grid topology as well as the power injections $P$ $P$ and thus change every hour. Hence, the weak spots of a grid also vary.

The situation becomes more involved if we go beyond the linear dc approximation. For instance, voltage instabilities played an important role during the blackout in the western USA in July 1996, where strongly nonlocal impacts were observed (Venkatasubramanian and Li, 2004; Hines, Dobson, and Rezaei, 2017).

3. Influence graphs

Influence graphs were proposed by Hines, Dobson, and Rezaei (2017) to describe the propagation of cascades and to understand nonlocal effects. This description was based on the statistical analysis of a large number of cascading failures simulated with a detailed physical model of the grid. The simulation data are used to determine the conditional probability $R j, i, m$ $R_{j, i, m}$ that a line $j$ $j$ fails in generation $m + 1$ $m + 1$ of a cascade given that line $i$ $i$ failed in generation $m$ $m$ . It turns out that probabilities do not change much with $m$ $m$ as long as $m \geq 1$ $m \geq 1$ . Only the first generation $m = 0$ $m = 0$ is different, as the grid in the initial state is typically $N - 1$ $N - 1$ secure such that failures are less probable. Hence, cascading failures are readily characterized using two matrices $R 0, R 1 \in R L \times L$ $R_{0}, R_{1} \in R^{L \times L}$ with elements $R j, i, 0$ $R_{j, i, 0}$ and $R j, i, 1 +$ $R_{j, i, 1 +}$ , respectively. An elementary example of a power grid and the reconstructed influence graph is shown in Fig. 24. A refined treatment including multiple line outages as separate states in a Markov chain model was proposed by Zhou et al. (2020).

Influence graphs are used to describe the propagation of cascading failures and to identify critical infrastructures. Left panel: power flows in an elementary test grid consisting of six nodes. Right panel: influence graph summarizing the likelihood that the failure of one line will induce the failure of another line. The nodes of this influence graph correspond to the lines in the original grid as identified by their end points. Adapted from Hines, Dobson, and Rezaei, 2017.

The influence graph can be used to identify critical lines as follows. If $P i, m$ $P_{i, m}$ is the probability that line $i$ $i$ fails in generation $m$ $m$ , then in the following generation

P j, m + 1 = \sum i R j, i, m P i, m .

$P_{j, m + 1} = \sum_{i} R_{j, i, m} P_{i, m} .$

(153)

Upon the introduction of a vectorial notation $P m = (P 1, m, \dots, P L, m) ⊤$ $P_{m} = (P_{1, m}, \dots, P_{L, m})^{⊤}$ , the probabilities evolve as $P 1 = R 0 P 0$ $P_{1} = R_{0} P_{0}$ and $P m + 1 = R 1 P m$ $P_{m + 1} = R_{1} P_{m}$ for all the following steps $m = 1, 2, \dots$ $m = 1, 2, \dots$ . The probability that a line $j$ $j$ will fail at some time during a long cascade is thus given by

P \infty = P 0 + \infty \sum m = 0 R m 1 R 0 P 0 .

$P_{\infty} = P_{0} + \sum_{m = 0}^{\infty} R_{1}^{m} R_{0} P_{0} .$

We can now analyze which elements have the largest impact on cascading failures. To this end one can calculate how the cumulative failure rate $P tot = \sum j P \infty, j$ $P_{tot} = \sum_{j} P_{\infty, j}$ will change if a single row $i$ $i$ of the matrices is $R 0$ $R_{0}$ and $R 1$ $R_{1}$ is modified, representing an upgrade of the respective line $i$ $i$ . One finds that most upgrades have little effect, but some lines lead to a drastic change of $P tot$ $P_{tot}$ . The corresponding lines represent the critical lines of the grid.

4. Statistical analysis of cascading failures

Sections VII.C.1–VII.C.3 have provided a microscopic picture of how flows are rerouted and how cascades propagate in electric power grids. But how do these aspects manifest in real-world large-scale grids under various operating conditions? A large-scale statistical analysis of these questions was recently provided by Yang and Motter (2017) using a high quality model of the three North American synchronous power grids (Eastern, Western, and Texas Interconnection). The utilized cascade model was significantly extended in comparison to the elementary model of Sec. VII.C.1, including load shedding, generator adaptation, and a model for transmission line overheating. A large number of cascades were simulated for a variety of different snapshots of the grid: different points in time with generation and load patterns. All cascades were triggered by the simultaneous failure of three lines randomly selected from the entire grid (Western, Texas) or from a certain grid area (Eastern Interconnection). Yang and Motter (2017) then evaluated the probability $P (p) ℓ$ $P_{ℓ}^{(p)}$ that a line $ℓ$ $ℓ$ will fail during a cascade as well as the probability $P (s) ℓ$ $P_{ℓ}^{(s)}$ that a line $ℓ$ $ℓ$ will be disconnected from the grid while carrying no load after the cascade.

The first important result of the statistical analysis is that coreness is an important factor that determines the probabilities $P (p, s) ℓ$ $P_{ℓ}^{(p, s)}$ ; see Fig. 25. Links with a higher coreness are more strongly connected and thus offer more options for flow rerouting. As a consequence they are more likely to suffer failures (the fraction of links with $P (p) ℓ > 0$ $P_{ℓ}^{(p)} > 0$ increases with the coreness) and they also have a higher average value of $⟨ P (p) ℓ ⟩$ $⟨ P_{ℓ}^{(p)} ⟩$ . In contrast, the probability of becoming disconnected (i.e., the fraction of links with $P (s) ℓ > 0$ $P_{ℓ}^{(s)} > 0$ ) decreases with the coreness. Links with coreness 1 can be disconnected by a single failure, while links with higher coreness have multiple connections to the grid and are thus hardly disconnected. Notably the average $⟨ P (s) ℓ ⟩$ $⟨ P_{ℓ}^{(s)} ⟩$ still increases with the coreness. In this analysis it should be noted that the vast majority of all links have a coreness of 2.

Statistical analysis of cascading failures. The coreness of a link determines its vulnerability to fail (upper panels) or to be disconnected (lower panels) in a cascade of failures. Shown is the fraction of lines that (c) fail or (e) get disconnected in at least one cascade. (d) Probability of failure $⟨ P (p) ℓ ⟩$ $⟨ P_{ℓ}^{(p)} ⟩$ . (f) Probability of disconnection $⟨ P (c) ℓ ⟩$ $⟨ P_{ℓ}^{(c)} ⟩$ averaged over all cascades and operating conditions. Simulations were carried out for the three North American grids under various operating conditions; see the text for details. Adapted from Yang and Motter, 2017.

In a second step, Yang and Motter (2017) investigated the vulnerable lines of the grids, i.e., lines with a high probability of failure ( $P (p) ℓ > 0.0005$ $P_{ℓ}^{(p)} > 0.0005$ ). It was found that the set of vulnerable lines is surprisingly small, including, for instance, only 48 out of 7637 lines in the Texas Interconnection. Moreover, there is a significant overlap of the vulnerable sets between different grid snapshots. That is, generation and load patterns determine the vulnerability only to a limited extent: many lines are vulnerable for various operating conditions. These findings raise the hope that grid stability can be significantly extended with comparably few grid extensions.

Finally, we are led to the question of which trigger events are particularly critical to grid stability. To answer this question Yang and Motter (2017) analyzed all cascades causing a load shedding above 300 MW. It was found that in these cases the three trigger events are typically close to each other and close to the vulnerable set of the grid, in terms of both geodesic network distances and geographic distances. This corresponds to our prior findings that flow rerouting is predominantly local.

5. The size of cascading failures

Cascades of failures can cause large-scale power outages. Statistics of the size of power outages are available for a variety of power grids, as shown in Fig. 26 for the North American power system. One observes that large-scale outages are actually not rare events; that is, the likelihood $P$ $P$ of an outage vanishes slowly with the size of the outage $N out$ $N_{out}$ . Indeed, it has been claimed that the distribution shows a power-law behavior

P \sim N - β out

$P \sim N_{out}^{- β}$

for large $N out$ $N_{out}$ with an exponent $1.3 \leq β \leq 2.0$ $1.3 \leq β \leq 2.0$ for various grids (Dobson et al., 2007). While it is generally hard to strictly establish a power law from empirical data (Clauset, Shalizi, and Newman, 2009), the statistics indicate that large-scale outages occur regularly.

Statistics of power outages in North American power grids (1984–1998). The empirical probability of a power outage as a function of the outage size on a logarithmic scale. The data suggest a power-law distribution with exponent $1.3 \leq β \leq 2.0$ $1.3 \leq β \leq 2.0$ . From Dobson *et al.*, 2007.

A detailed analysis of cascading failures including empirical statistics was provided by Dobson et al. (2007). They analyzed outage statistics both in detailed numerical models and in a coarse-scaled model admitting a closed form solution (Dobson, Carreras, and Newman, 2005). This model reduces a cascade to the following three essential elements:

• Many elements with heterogeneous loading.—The model considers $N$ $N$ infrastructure elements whose initial load $L 0$ $L_{0}$ is drawn at random from a uniform distribution in the interval $[0, L max]$ $[0, L_{\max}]$ . For simplicity we normalize all quantities such that $L max = 1$ $L_{\max} = 1$ .

• Cascade mechanism.—If the loading of an element exceeds a limit $L \geq L fail$ $L \geq L_{fail}$ , then this element fails and the loading is distributed to other elements. For simplicity we assume that the loading of all other elements then increases by the value $Δ L 1$ $Δ L_{1}$ .

• Initial trigger.—The cascade is triggered by some initial failures, thus increasing the initial load of all elements by an amount $Δ L 0$ $Δ L_{0}$ .

For $L fail = 1$ $L_{fail} = 1$ , $Δ L 0 + Δ L 1 N \leq 1$ $Δ L_{0} + Δ L_{1} N \leq 1$ , the number of failing elements $N out$ $N_{out}$ follows a quasibinomial distribution

P (N out = k) = (N k) Δ L 0 (Δ L 0 + k Δ L 1) k - 1 \times (1 - Δ L 0 - k Δ L 1) N - k .

$P (N_{out} = k) = (\binom{N}{k}) Δ L_{0} (Δ L_{0} + k Δ L_{1})^{k - 1} \times (1 - Δ L_{0} - k Δ L_{1})^{N - k} .$

(154)

If $Δ L 1 ≪ 1 / N$ $Δ L_{1} ≪ 1 / N$ , this distribution strongly peaks around the mean. But if $Δ L 1$ $Δ L_{1}$ increases and approaches the value $1 / N$ $1 / N$ , the system approaches a critical point with a power-law distribution

P (N out = k) \sim k - 1.5,

$P (N_{out} = k) \sim k^{- 1.5},$

(155)

thereby providing a fair fit to the empirical results.

But why should a power system remain critical over many years? Dobson et al. (2007) argued that criticality emerges in a self-organized way. The loads in the system are increasing year to year; in terms of the model $L fail$ $L_{fail}$ decreases and $Δ L 1$ $Δ L_{1}$ increases. As the system approaches the critical point, the probability and size of power outages increases and countermeasures are implemented, such as via an extension of infrastructures. Hence, the system remains close to the critical point. This behavior was reproduced in a detailed numerical model (Carreras et al., 2002) in which power dispatch and cascades occur on short timescales (cf. Secs. III.B and VII.C.1) and load and infrastructures evolve on much slower timescales. In the model, the slow dynamics brings the power system close to criticality and such large-scale cascades become likely. Note that the hypothesis of self-organized criticality is highly controversial, as summarized by Fairley (2004). Different explanations were put forward by Nesti, Sloothaak, and Zwart (2020) relating power-law distributions in outages sizes to power-law distributions in city sizes, and by Kaiser and Witthaut (2021b).

6. Transient effects in cascading failures

Our previous analysis of cascading failures was restricted to a quasistatic picture. We considered the existence and the properties of the steady state after the failure of a transmission line while neglecting any dynamic. However, the existence of a stable state after a line failure is only a necessary condition for grid stability, not a sufficient one. It is not a priori clear that the grid will end up in this state, but there are also other possibilities: (i) It is not guaranteed that the system relaxed to a fixed point after failure. Instead, it can also relax to a limit cycle or another attractor, as discussed in Sec. V.C.3. (ii) Even if the system relaxes to the desired steady state, it may violate operational constraints during the relaxation, which might also cause secondary outages. This relates to the aspect of survivability introduced in Sec. V.C.3. As a consequence, transient effects can increase the vulnerability of grid and a quasistatic picture can miss important threats to stability (Simonsen et al., 2008).

A detailed analysis of cascading failures, including transient overloads, was provided by Schäfer, Witthaut et al. (2018). The dynamics was modeled in terms of the aggregated dynamical model of Sec. III.C.9 and a line was instantaneously removed from the simulation when the flow exceeded a line rating $| F i j (t) | > F max i j$ $| F_{i j} (t) | > F_{i j}^{\max}$ . The impact of such a transient outage is illustrated for an elementary test grid in Fig. 27, where a cascade is triggered by the failure of a single line $(2, 4)$ $(2, 4)$ [Fig. 27(a)]. In a quasistatic picture, the grid immediately jumps to a new stable steady state [Fig. 27(c)]. However, the line ratings are exceeded transiently, leading to secondary failures of line $(4, 5)$ $(4, 5)$ and eventually to a cascade disconnecting the grid.

Transient effects can increase the vulnerability of a grid to cascading failures. (a) An elementary test grid with two effective generators ( $P = + 1.5 s - 2$ $P = + 1.5 s^{- 2}$ ; squares) and three effective consumers ( $P = - 1 s - 2$ $P = - 1 s^{- 2}$ ; circles). All lines have $K i j = 10 / 6$ $K_{i j} = 10 / 6$ and $F max i j = 1$ $F_{i j}^{\max} = 1$ . At time $t = 1 s$ $t = 1 s$ , the line $(2, 4)$ $(2, 4)$ fails. (b) In a full dynamical study, line ratings are exceeded transiently, leading to secondary line outages. Shown is the number of failed lines as a a function of time. (c) In a quasistatic picture, the grid immediately relaxes back to a stable steady state that respects the line ratings. Shown are the steady-state line flows $| F ℓ |$ $| F_{ℓ} |$ before and after the initial failure. (d) Time evolution of the line flows $| F ℓ (t) |$ $| F_{ℓ} (t) |$ taking into account both the continuous dynamics and the outage of lines. From Schäfer, Witthaut *et al.*, 2018.

Whether transient overloads are important depends crucially on grid properties, the line ratings $F max i j$ $F_{i j}^{\max}$ , and the initial trigger that starts the cascade. The strongest impacts have been observed for heavily loaded grids, where the likelihood of secondary overloads more than doubles.

D. Braess’s paradox

The loss of a transmission line can cause a blackout, either directly, as discussed in Sec. V.A, or indirectly via a cascade of failures. But the reinforcement of a transmission line or the addition of a new line can also induce a loss of stability (Witthaut and Timme, 2012; Coletta and Jacquod, 2016). This surprising effect is demonstrated for an elementary test grid in Fig. 28. The initial grid relaxes to a synchronized fixed point, whereas synchronization becomes impossible after certain grid extensions and the synchronized fixed point has ceased to exist.

Braess’s paradox: loss of the synchronous fixed point due to grid extensions. (a)–(c) Topology of the network. The vertices generate or consume the power $P i = \pm P$ $P_{i} = \pm P$ , and the transmission lines have a capacity $K$ $K$ . We consider the cases where (b) the capacity of one line is doubles and (c) a new line is added. (d)–(f) Dynamics of the grid given by the equations of motion (111) starting from the initial state $θ i (0) = ˙ θ i (0) = 0$ $θ_{i} (0) = {\dot{θ}}_{i} (0) = 0$ . The initial grid relaxes to a synchronized fixed point, while synchronization becomes impossible after the grid extensions. The parameters are $K = 1.03 P$ $K = 1.03 P$ , $D = P$ $D = P$ , and $J = 1$ $J = 1$ . Adapted from Witthaut and Timme, 2012.

We can understand the loss of a fixed point using the cycle flow approach introduced in Sec. VI.C.3. Increasing the capacity of a line as in Fig. 28(b) will always lower the load on this specific line and thus relieve the grid. However, any fixed point must also satisfy the winding number constraint (124) for every cycle in the network. As a consequence, all flows in the grid $F$ $F$ are affected, which may cause the disappearance of synchronized phase-locked states, a prerequisite for normal grid operation. In particular, for the example shown in Fig. 28(b) a cycle flow must be added in the counterclockwise direction to continue to satisfy the constraint (124) after the line extension. But this would increase the load on lines $4 \to 5$ $4 \to 5$ and $4 \to 8$ $4 \to 8$ above the upper limit $K$ $K$ , which is not possible. The fixed point ceases to exist.

Braess’s paradox can be used to improve grid stability. In certain situations it can be beneficial to shut down a transmission line on purpose to relieve another line (Motter, 2004; Witthaut and Timme, 2015); cf. Sec. VII.A.5. Similar effects were observed for other supply networks, starting with a seminal work on traffic networks by Braess (1968).

VIII. Perturbations, Fluctuations, and Transient Dynamics

Since the ascent of renewable energy sources, accompanied by increased trading and regulatory activities, power inputs and outputs increasingly fluctuate on all timescales. In addition, consumption patterns keep changing due to the increasing use of smart devices, digitalization, and global connectivity. Such perturbations and fluctuations require analyses of the state of a power grid as a driven, nonequilibrium system where voltages, voltage angles, and flows are nonstatic and even nonstationary as they respond dynamically to time-varying signals.

In this section, we address key questions about perturbations and fluctuations in power grids. In a mathematical model for bulk grids, we illustrate in Sec. VIII.A how fluctuations of feed-in and consumed power translate to frequency fluctuations while highlighting their non-Gaussian statistics. We also highlight the impact of grid heterogeneities on non-Gaussian features and, in Sec. VIII.B, present driving response relations for networked systems that explicitly take into account any given grid topology. Beyond deriving the general form of linear response theory, we explain collectively emerging dynamic phenomena including resonances, bulk oscillations, localization, and signal propagation in Sec. VIII.C. In Sec. VIII.D, we outline how the Wentzel-Kramers-Brillouin (WKB) method helps one to predict the probability of blackouts if rare but large fluctuations kick the system out of a stable phase-locked state. For fluctuations in production and consumption occurring simultaneously all over the grid, a natural implementation in the spirit of statistical mechanics amounts to considering ensembles of power grids (Sec. VIII.E) with possible applications to obtain macroscopic average quantities such as the market volume and market costs.

A. Frequency fluctuations from time series data

Dynamic driving signals continually change the overall state of a power grid. Because of fluctuations, the collective dynamics of electric power grids thereby is not only intrinsically out of equilibrium but also intrinsically nonstationary. Fluctuations often directly modify essential parameters of input and output and make them time dependent, reflecting fluctuating power feed-in and consumption. One key example is the fluctuating power occurring as a time-dependent parameter $P i (t)$ $P_{i} (t)$ , such as in the second-order model (111).

To analyze how the statistics of input power fluctuations influence the grid frequency, we first focus on the bulk frequency dynamics given by Eq. (87), assuming that the power imbalance is due to the following stochastic fluctuations at the grid nodes $i = 1, \dots, N$ $i = 1, \dots, N$ with strength $q i$ $q_{i}$ :

Δ P = N \sum i = 1 q i ξ i (t) ≕ ¯ q ¯ ξ (t) .

$Δ P = \sum_{i = 1}^{N} q_{i} ξ_{i} (t) ≕ \bar{q} \bar{ξ} (t) .$

(156)

In traditional analyses of power engineering, input fluctuations are either neglected or modeled to be of a Gaussian nature, with either white or colored noise, i.e., without any or with small temporal correlations; see Zhang and Li (2010), Wood, Wollenberg, and Sheblé (2014), and Schäfer et al. (2017). For Gaussian white noise power fluctuations, the resulting distribution of frequency deviations $¯ ω$ $\bar{ω}$ is Gaussian with zero average and a standard deviation $σ = \sqrt \sum i q i / 2 η ¯ J 2$ $σ = \sqrt{\sum_{i} q_{i} / 2 η {\bar{J}}^{2}}$ . Such a Gaussian model neglects heavy tails observed in bulk frequency distributions (Schäfer, Beck et al., 2018; Gorjão, Jumar et al., 2020). As Fig. 29 illustrates, such heavy tails indicate that strong deviations from the reference frequency are orders of magnitude more frequent than a Gaussian model would predict.

Non-Gaussian frequency deviations with a reduced number of extreme events. (a) Distribution of frequencies around a set frequency (50 Hz) in the Continental European grid (yellow) deviates from the best Gaussian distribution (parabolic blue solid line). In particular, the probability (per Hz) of large events of more than 0.1 Hz frequency deviations is increased by more than a factor of 10 compared to the Gaussian fit. (b) Counts of “extreme events,” i.e., the number of frequency deviations larger than 0.1 Hz resolved for each minute, within full hours, accumulated for 2017 and 2011, respectively. Compiled using data from Reseau de Transport d’Electricite (2017) and the methods of Schäfer, Timme, and Witthaut (2018).

Heavy-tailed frequency distributions may emerge from either heavy-tailed distributions of the power fluctuations $¯ ξ (t)$ $\bar{ξ} (t)$ or the temporal variability of the system parameters (Schäfer, Beck et al., 2018). In particular, abrupt changes of power generation at the beginning of the trading intervals may cause large frequency deviations (Gorjão, Anvari et al., 2020). Counting the number of extreme events of the frequencies, defined as frequency deviations of more than 100 mHz from the set frequency, and comparing these data from the years 2017 and 2011 reveal additional information [Fig. 29(b)]; cf. Schäfer, Timme, and Witthaut (2018). First, the total number of extreme events is reduced in 2017. Second, the number of these events is predominantly reduced in an interval of a few minutes after full hours. Third, the total number of threshold violations seems to still be substantial in the first few minutes, as well as in about the 18th, the 33rd, and the 48th minutes after a full hour (each a few minutes after a full quarter hour), hinting that switching the trading intervals from each hour (in 2011) to each quarter of an hour (in 2017), and thereby reducing the volume of trading per event together with changes in regulatory action, might have caused the reduction in the number of extreme frequency deviations. The histogram also hints at a characteristic timescale of response of frequency deviations at approximately 2–4 min delay. Both the details of such deviations and the generality of their occurrence remain poorly understood to date.

The resulting extreme events pose theoretical questions for analysis and serious practical challenges, such as for security assessment. Recent modeling work (Wolff et al., 2019) illustrated that in power grids in which consumption, generation, and transmission infrastructures are heterogeneous, fluctuating wind power injection at nodes that are weakly coupled to the grid particularly contribute non-Gaussian features to frequency deviations. Some further studies considered non-Gaussian effects, focusing on either theoretical aspects of how they may emerge or numerical evaluations of individual wind and solar data; see Kashima, Aoyama, and Ohta (2015), Anvari et al. (2016), Schmietendorf, Peinke, and Kamps (2017), and Totz, Olmi, and Schöll (2020). In particular, the intermittent nature of short-term wind fluctuations cause novel types of frequency and voltage fluctuations and thereby influence stability properties of grid dynamics (Schmietendorf, Peinke, and Kamps, 2017). Wind power induced fluctuations, quantified in terms of frequency deviation variance, moreover, propagate along interconnected synchronous machines in a characteristic way, with exponentially decaying amplitudes (Haehne et al., 2019).

B. Network linear response theory for fluctuating power

How does the collective grid dynamics respond to input fluctuations? We review the general linear response theory valid for small perturbations with arbitrary time dependence (Haehne et al., 2019; Zhang et al., 2019; Zhang, Witthaut, and Timme, 2020). Consider a small perturbation $δ P (t) = [δ P 1 (t), \dots, δ P N (t)] T$ $δ P (t) = [δ P_{1} (t), \dots, δ P_{N} (t)]^{T}$ that at all nodes $i$ $i$ results in phase deviations $α i (t)$ $α_{i} (t)$ To first order in the deviation $α i$ $α_{i}$ , perturbations evolve according to the linear wave equation (112), as governed by the weighted graph Laplace matrix $L \in R N \times N$ $L \in R^{N \times N}$ , which is defined in Eq. (113). For example, for homogeneous parameters the inertia $J i = J$ $J_{i} = J$ , and the damping factor $D i = D$ $D_{i} = D$ , the response of the phase $α i (t)$ $α_{i} (t)$ to dynamic perturbations, such as a change of power $δ P j (t)$ $δ P_{j} (t)$ at node $j$ $j$ , is given by (Haehne et al., 2019)

α i (t) = J D 2 ω R \int t - \infty \sum j δ P j (t') G i j (t' - t) d t' τ,

$α_{i} (t) = \frac{J}{D^{2} ω_{R}} \int_{- \infty}^{t} \sum_{j} δ P_{j} (t^{'}) G_{i j} (t^{'} - t) \frac{d t^{'}}{τ},$

(157)

with the grid reference frequency $ω R$ $ω_{R}$ ; see Sec. III.C.1. The propagator from node $i$ $i$ to node $j$ $j$ is defined by

G i j (t' - t) = N \sum n = 1 \sum σ = \pm 1 ϕ n i ϕ * n j 2 \sqrt 1 - ˜ Λ n (- σ) e (1 + σ \sqrt 1 - ˜ Λ n) (t' - t) / τ .

$G_{i j} (t^{'} - t) = \sum_{n = 1}^{N} \sum_{σ = \pm 1} \frac{ϕ_{n i} ϕ_{n j}^{*}}{2 \sqrt{1 - {\tilde{Λ}}_{n}}} (- σ) e^{(1 + σ \sqrt{1 - {\tilde{Λ}}_{n}}) (t^{'} - t) / τ} .$

(158)

See Haehne et al. (2019) for the full derivation. In Eq. (158) $1 / τ = D / J$ $1 / τ = D / J$ is the local relaxation rate (of a single node disconnected from other nodes) and $˜ Λ n = (J / D 2 ω R) Λ n$ ${\tilde{Λ}}_{n} = (J / D^{2} ω_{R}) Λ_{n}$ , where $Λ n \in R$ $Λ_{n} \in R$ are eigenvalues and $ϕ n \in C N$ $ϕ_{n} \in C^{N}$ are the corresponding eigenvectors of the generalized graph Laplacian matrix $L$ $L$ [Eq. (113)], which is related to the stability matrix used in small signal stability analysis (Milano, 2010; Zhang, Rehtanz, and Pal, 2012). For each eigenvalue of the generalized graph Laplacian $Λ n \in R$ $Λ_{n} \in R$ there are two eigenvalues of the linearized swing equations (84), the linear wave equation (112), as given by $ϵ n σ = - i (1 + σ \sqrt 1 - ˜ Λ n) / τ$ $ϵ_{n σ} = - i (1 + σ \sqrt{1 - {\tilde{Λ}}_{n}}) / τ$ , where $σ = \pm$ $σ = \pm$ . Thus, the denominator in Eq. (158) is proportional to the difference of these eigenvalues $2 \sqrt 1 - ˜ Λ n = τ i (ϵ n + - ϵ n -)$ $2 \sqrt{1 - {\tilde{Λ}}_{n}} = τ i (ϵ_{n +} - ϵ_{n -})$ . Note that the relaxation rate of each mode is given by $Γ n σ = - Im ϵ n σ$ $Γ_{n σ} = - Im ϵ_{n σ}$ , which is identical to the local relaxation rate $1 / τ$ $1 / τ$ if $˜ Λ n > 1$ ${\tilde{Λ}}_{n} > 1$ but differs from it by the term $σ \sqrt 1 - ˜ Λ n / τ$ $σ \sqrt{1 - {\tilde{Λ}}_{n}} / τ$ if $˜ Λ n < 1$ ${\tilde{Λ}}_{n} < 1$ , yielding for the mode $σ = - 1$ $σ = - 1$ a slower relaxation in the grid than for individual nodes. This expression applies to any disturbance $δ P j (t)$ $δ P_{j} (t)$ as long as its amplitude is sufficiently small; see Eq. (S2) in the Supplemental Material of Tamrakar, Conrath, and Kettemann (2018) for a validity condition. Similar expressions can be obtained in response to a change of any of the system parameters, like a change of power capacitance $δ K i j (t)$ $δ K_{i j} (t)$ between the nodes $i$ $i$ and $j$ $j$ (Kettemann, 2016; Witthaut et al., 2016; Manik et al., 2017).

In the frequency representation of the linear response, we write the phase deviation $α i (t)$ $α_{i} (t)$ as a generalized Fourier series and expand its spatial dependence in terms of eigenvectors $ϕ n$ $ϕ_{n}$ of the generalized Laplacian $L$ $L$ . Thereby we obtain (Kettemann, 2016; Auer et al., 2017; Tamrakar, Conrath, and Kettemann, 2018; Zhang et al., 2019)

α i (t) = \int \infty - \infty N \sum n = 1 c n (ε) ϕ n i e - i ε t d ε,

$α_{i} (t) = \int_{- \infty}^{\infty} \sum_{n = 1}^{N} c_{n} (ε) ϕ_{n i} e^{- i ε t} d ε,$

(159)

where $c n (ε)$ $c_{n} (ε)$ is the contribution strength of the angular frequency $ε$ $ε$ at the $n$ $n$ th node. Likewise expanding the disturbance in a Fourier series, we obtain

δ P i (t) = D 2 ω R J \int \infty - \infty N \sum n = 1 η n (ε) ϕ n i e - i ε t d ε .

$δ P_{i} (t) = \frac{D^{2} ω_{R}}{J} \int_{- \infty}^{\infty} \sum_{n = 1}^{N} η_{n} (ε) ϕ_{n i} e^{- i ε t} d ε .$

(160)

Inserting these expansions for $α i (t)$ $α_{i} (t)$ and $δ P i (t)$ $δ P_{i} (t)$ into the linear wave equation (112), one finds, requiring that the equation is fulfilled for each term of the Fourier series, $(- τ 2 ε 2 - i 2 τ ε + ˜ Λ n) c n (ε) = η n (ε)$ $(- τ^{2} ε^{2} - i 2 τ ε + {\tilde{Λ}}_{n}) c_{n} (ε) = η_{n} (ε)$ . For a given disturbance, the Fourier component of the phase deviation $c n (ε)$ $c_{n} (ε)$ is thus given in response to that of the disturbance $η n (ε)$ $η_{n} (ε)$ . Inserting that expression for $c n (ε)$ $c_{n} (ε)$ back into the Fourier series, one gets

α i (t) = \int \infty - \infty N \sum n = 1 (- τ 2 ε 2 - i 2 τ ε + ˜ Λ n) - 1 η n (ε) ϕ n i e - i ε t d ε .

$α_{i} (t) = \int_{- \infty}^{\infty} \sum_{n = 1}^{N} {(- τ^{2} ε^{2} - i 2 τ ε + {\tilde{Λ}}_{n})}^{- 1} η_{n} (ε) ϕ_{n i} e^{- i ε t} d ε .$

(161)

A linear response theory generalizing Eq. (157) to include inhomogeneous parameters and the presence of Ohmic losses was derived by Plietzsch et al. (2019), and the case without inertia ( $J = 0$ $J = 0$ ) was analyzed by Tyloo, Coletta, and Jacquod (2018). Zhang et al. (2019) disentangled the responses and analyzed spatiotemporal response patterns (to be discussed in Sec. VIII.C) starting with a focus of driving one node at one frequency, i.e., $η n (ε) = δ (ε - ε 0) δ n, j$ $η_{n} (ε) = δ (ε - ε_{0}) δ_{n, j}$ in Eq. (160), to obtain characteristics of linear response estimates analytically.

Given the phase deviation $α i (t)$ $α_{i} (t)$ for the time $t$ $t$ and the position $i$ $i$ , one can calculate the temporal and spatial evolution of phase deviations, the frequency deviations $(δ f) i (t) = \partial t α i (t)$ $(δ f)_{i} (t) = \partial_{t} α_{i} (t)$ , and the rate of change of the frequency deviations $\partial t (δ f i (t))$ $\partial_{t} (δ f_{i} (t))$ (Pagnier and Jacquod, 2019a), as well as time-averaged moments of these quantities. Furthermore, moments of increments of frequency at the node $i$ $i$ , $f i Δ t = f i (t + Δ t) - f i (t)$ $f_{i Δ t} = f_{i} (t + Δ t) - f_{i} (t)$ (Haehne et al., 2019), contain information about correlations at timescale $Δ t$ $Δ t$ . In the following, we review the results obtained in linear response for these quantities for different types of nonstationary signals impinging on exemplary power grid structures.

C. Spatiotemporal responses from localized to resonant

1. Propagation of short-duration disturbances

For disturbances that last for a short time compared to the local relaxation time $τ = J / D$ $τ = J / D$ , the linear response at the node $i$ $i$ due to a power pulse

δ P k (t) = δ P δ k j τ δ (t - t 0)

$δ P_{k} (t) = δ P δ_{k j} τ δ (t - t_{0})$

(162)

at the time $t 0$ $t_{0}$ and the node $j$ $j$ by Eq. (157) yields, for the times $t > t 0$ $t > t_{0}$ ,

α i (t) = - J δ P 2 D 2 ω R N \sum n = 1 \sum σ = \pm 1 ϕ n i ϕ * n j \sqrt 1 - ˜ Λ n e - (1 + σ \sqrt 1 - ˜ Λ n) (t - t 0) / τ .

$α_{i} (t) = - \frac{J δ P}{2 D^{2} ω_{R}} \sum_{n = 1}^{N} \sum_{σ = \pm 1} \frac{ϕ_{n i} ϕ_{n j}^{*}}{\sqrt{1 - {\tilde{Λ}}_{n}}} e^{- (1 + σ \sqrt{1 - {\tilde{Λ}}_{n}}) (t - t_{0}) / τ} .$

(163)

For meshed, spatially embedded grids and sufficiently large inertia, analytical (Kettemann, 2016) and numerical results (Tamrakar, Conrath, and Kettemann, 2018) showed that disturbances propagate ballistically with the velocity $v$ $v$ such that the time $t > t 0$ $t > t_{0}$ when signals arrive at the geometrical distance $r$ $r$ from the position of the disturbance is given by

t - t 0 = r / v .

$t - t_{0} = r / v .$

(164)

In Eq. (164) the arrival time $t > t 0$ $t > t_{0}$ is defined as the time that it takes after a single-node disturbance in $δ P$ $δ P$ at time $t 0$ $t_{0}$ until a phase deviation $α i (t)$ $α_{i} (t)$ exceeds a threshold $α th$ $α_{th}$ (in the example in Fig. 30, $α th = 10 - 6 J δ P / D 2 ω R$ $α_{th} = 10^{- 6} J δ P / D^{2} ω_{R}$ is chosen).

Arrival time $t * = (t - t 0) / τ$ $t^{*} = (t - t_{0}) / τ$ of a short pulse disturbance in units of $τ = J / D$ $τ = J / D$ vs the geometric distance $r * = r / a$ $r^{*} = r / a$ in units of the transmission line length $a$ $a$ . The green rhombi represent the numerical results for (a) a square grid and (b) German transmission grid topology: $K J / D 2 ω R = 105$ $K J / D^{2} ω_{R} = 10^{5}$ , a realistic value for high voltage transmission grids. Also shown are the ballistic equation (164) with the fitted velocity $v$ $v$ (red line) and the analytic velocity equation (165) (red dashed line). The black squares indicate node-resolved arrival times $t * i$ $t_{i}^{*}$ . Unfilled squares indicate arrival times, averaged over all nodes at the same distance (c) for a Cayley tree with branching $b = 2$ $b = 2$ , and (d) for the German transmission grid. In (c) and (d) $K J / D 2 ω R = 10$ $K J / D^{2} ω_{R} = 10$ . Diffusive equation (166) fitted (pink line) and with the analytical equation (167) (pink dashed line). Also shown are the ballistic equation (164) fitted (red line) and with analytic velocity (red dashed line). $σ = P / P c = 0.1$ $σ = P / P_{c} = 0.1$ , with critical power $P c$ $P_{c}$ for the respective grids. Adapted from Tamrakar, Conrath, and Kettemann, 2018.

For homogeneous regular grids, the velocity $v$ $v$ can be derived as follows from the response theory equation (163) (Kettemann, 2016; Tamrakar, Conrath, and Kettemann, 2018):

v = a \sqrt K J ω R (1 - P 2 P 2 c) 1 / 4,

$v = a \sqrt{\frac{K}{J ω_{R}}} {(1 - \frac{P^{2}}{P_{c}^{2}})}^{1 / 4},$

(165)

with the length of a single transmission line $a$ $a$ . Here, $P < P c$ $P < P_{c}$ , where $P c$ $P_{c}$ is the critical power above which no stationary solution to the power flow equations exists. Equation (165) gives lower bounds for the arrival times [Eq. (164)] for both regular square and German transmission grid topology [the red dashed lines in Figs. 30(a) and 30(b)]. According to Eq. (165) the maximal velocity of disturbances in both regular and meshed real grid topologies increases with the decreasing inertia $J$ $J$ and decrease as $P \to P c$ $P \to P_{c}$ .

For unmeshed grid topologies it was found by Tamrakar, Conrath, and Kettemann (2018) that arrival times grow not linearly but quadratically with distance $r$ $r$ , resembling the diffusion

t - t 0 = r 2 / 4 D .

$t - t_{0} = r^{2} / 4 D .$

(166)

For Cayley tree graphs with branching number $b$ $b$ the diffusion constant is derived as

D (b) = τ Δ 2 \sqrt b \sqrt b - 1,

$D (b) = \frac{τ Δ^{2} \sqrt{b}}{\sqrt{b} - 1},$

(167)

where $Δ = (K / J ω R) 1 / 2 (1 - P 2 / P 2 c) 1 / 4 (\sqrt b - 1)$ $Δ = (K / J ω_{R})^{1 / 2} (1 - P^{2} / P_{c}^{2})^{1 / 4} (\sqrt{b} - 1)$ is the positive Fiedler value, the spectral gap of the graph Laplacian, independent of the number of grid nodes $N$ $N$ . This is confirmed by the numerical calculations illustrated in Fig. 30(c) (dashed line). For low inertia ( $J < J c$ $J < J_{c}$ ) the collective dynamics of coupled nodes also results in diffusive spreading of the disturbances in meshed grids; see Fig. 30(d). $J c$ $J_{c}$ is obtained from the condition that slow modes with a small relaxation rate appear where the spectral gap (Fiedler value) $Δ$ $Δ$ is smaller than the local relaxation rate, ( $1 / τ = D / J$ $1 / τ = D / J$ ) (Tamrakar, Conrath, and Kettemann, 2018). As Fig. 30(d) illustrates, the spreading is more strongly delayed for some nodes, and disturbances are localized in certain grid regions where nodes do not become excited above the threshold (at least within the observation time). Linear response theory may explain this feature, as the response [Eq. (163)] is proportional to the eigenvector amplitude of the Laplacian $ϕ n i$ $ϕ_{n i}$ . As previously noted, the long-time transient behavior is dominated by the Fiedler vector, the eigenvector with the smallest nonzero eigenvalue, which is shown in Fig. 31 (left panel) for the German transmission grid topology and illustrates two geographic regions with high amplitudes. A similar result was obtained by Pagnier and Jacquod (2019a), who found the Fiedler vector of the weighted European transmission grid to be localized at the southern and northern borders. Moreover, the global RoCoF was found to decay with increasing system inertia, which they related to the Fiedler vector intensity. In random graphs, the Fiedler vector can be strongly localized, even on a single node [Fig. 31 (right panel)], so single-node disturbances of sufficiently small amplitudes remain localized. Thus, strong randomness and inhomogeneity may result in localization of disturbances, as noted by Kettemann (2016), a phenomenon well known as Anderson localization (Anderson, 1958).

Normalized Fiedler vector intensity $| ϕ 2, i | / max n | ϕ 2, n |$ $| ϕ_{2, i} | / \max_{n} | ϕ_{2, n} |$ relative to its maximum (shown in black) for the German transmission grid (left panel) featuring maximal intensity at the south and north borders and for a random grid with strong localization at a single node (black; at the bottom right of the right panel). Adapted Torres-Sánchez, Freitas de Abreu, and Kettemann, 2020.

2. Propagation of stochastic disturbances

On large timescales, frequency control measures compensate feed-in fluctuations of renewable generators, as reviewed in Sec. III.D, thereby maintaining stable grid operation. However, on timescales below 1 s, grid frequency fluctuations increase with increasing wind power production (Haehne et al., 2018). Moreover, the timescale separating local from interarea modes is also of the order of 1 s (Zhang, Rehtanz, and Pal, 2012). Are such fluctuations a local feature, for instance, resulting from locally high wind power injection, or do they affect grid dynamics over large ranges? To address this question, the subsecond grid frequency dynamics has been simulated by stochastically perturbing the grid. Model simulations of coupled nonlinear oscillator models with synthetically generated wind power feed-in time series (Haehne et al., 2019) indicate that the variance of short-term fluctuations decays for large inertia exponentially with distance to the feed-in node. These findings hold for both linear chain networks and German transmission grid topology (see Fig. 32), which is in agreement with analytical results for the variance of frequency increment distributions,

⟨ (^D Δ t ω i) 2 ⟩ = ⟨ (^D Δ t δ P 1) 2 ⟩ J K ω R exp (- d i, j ξ),

$⟨ ({\hat{D}}_{Δ t} ω_{i})^{2} ⟩ = \frac{⟨ ({\hat{D}}_{Δ t} δ P_{1})^{2} ⟩}{J K ω_{R}} \exp (- \frac{d_{i, j}}{ξ}),$

(168)

obtained by linear response theory [Eqs. (157) and (158)] for chainlike grids with $N ≫ 1$ $N ≫ 1$ nodes. Here $^D Δ t f (t) = f (t + Δ t) - f (t)$ ${\hat{D}}_{Δ t} f (t) = f (t + Δ t) - f (t)$ for any function $f$ $f$ , and $⟨ (^D Δ t δ P 1) 2 ⟩$ $⟨ ({\hat{D}}_{Δ t} δ P_{1})^{2} ⟩$ is the second moment of the increment distribution of the disturbance at site $j = 1$ $j = 1$ . Thus, the second moment of the frequency increments decays exponentially with topological distance $d i, j = i - 1$ $d_{i, j} = i - 1$ from the position of the disturbance with correlation length

ξ = v τ / 2 = \sqrt J K / 2 \sqrt ω R D .

$ξ = v τ / 2 = \sqrt{J K} / 2 \sqrt{ω_{R}} D .$

(169)

Below a critical inertia $J < J c = ω R D 2 N 2 / π 2 K \approx 1.6 \times 106 kgm 2$ $J < J_{c} = ω_{R} D^{2} N^{2} / π^{2} K \approx 1.6 \times 10^{6} {kgm}^{2}$ , there are nonzero eigenvalues $Λ n < 1$ $Λ_{n} < 1$ ; cf. Fig. 32. According to Eq. (158), the amplitude of these modes decays with a rate smaller than $1 / τ$ $1 / τ$ .

(a) Variance of increments $⟨^D Δ t ω 2 i ⟩$ $⟨ {\hat{D}}_{Δ t} ω_{i}^{2} ⟩$ for a chain of $N = 50$ $N = 50$ oscillators. The straight lines have a slope $m = - 1 / ξ$ $m = - 1 / ξ$ with Eq. (169). (b) Eigenvalues of Laplacian $Λ n$ $Λ_{n}$ for various inertia $J$ $J$ . (c) Fits of slope $m$ $m$ (orange) converge to an analytical $m = - 1 / ξ$ $m = - 1 / ξ$ [Eq. (169)] (blue) with increasing inertia $J$ $J$ . Vertical dashed line, $J c$ $J_{c}$ ; error bars, $2 σ$ $2 σ$ confidence bound. From Haehne *et al.*, 2019.

In sharp contrast, the kurtosis

k = ⟨ (^D Δ t f - ⟨^D Δ t f ⟩) 4 ⟩ ⟨ (^D Δ t f - ⟨^D Δ t f ⟩) 2 ⟩ 2

$k = \frac{⟨ ({\hat{D}}_{Δ t} f - ⟨ {\hat{D}}_{Δ t} f ⟩)^{4} ⟩}{⟨ ({\hat{D}}_{Δ t} f - ⟨ {\hat{D}}_{Δ t} f ⟩)^{2} ⟩^{2}}$

(170)

of frequency increments, quantifying deviations from Gaussian distribution ( $k = 3$ $k = 3$ ), is found to decay slowly, subexponentially, with distance from the disturbance. Thus, the non-Gaussian shape of frequency fluctuations (see Sec. VIII.A) persists over long ranges (Fig. 33).

Kurtosis $k$ $k$ of frequency increment PDF in a chain of $N = 100$ $N = 100$ nodes. The frequency increment distribution $p (^D Δ t ω i)$ $p ({\hat{D}}_{Δ t} ω_{i})$ deforms only slowly toward an almost Gaussian distribution. Insets (from left to right): $i \in {2, 20, 50}$ $i \in {2, 20, 50}$ , $θ = 0.01 s$ $θ = 0.01 s$ . From Haehne *et al.*, 2019.

In addition to these fundamental aspects, linear response arguments have been used in various applications. Using realistic grid models, analyses have addressed how fluctuations affect the primary control effort (Tyloo and Jacquod, 2021), where additional inertia should be placed (Pagnier and Jacquod, 2019b), and which fluctuation sources have the strongest impact on the grid (Gambuzza et al., 2017).

3. Localization, distributed resonances, and bulk oscillations

Which collective dynamical phenomena does linear response theory capture? Zhang et al. (2019) analyzed the response dependence on grid topology, on the exact location of perturbed and responding nodes in the network, and on the frequency content of the power fluctuations driving the system and compared the results to the full nonlinear system dynamics. If a given node $j$ $j$ is driven by power fluctuations at a given frequency $ω$ $ω$ , i.e., $η n (ε) = δ (ε - ω) δ n, j$ $η_{n} (ε) = δ (ε - ω) δ_{n, j}$ in Eq. (160), three regimes with qualitatively different stationary responses emerge; see Fig. 34.

Localization, distributed resonances, and bulk oscillations. In a small network example (a), the joint implications on response strength (b) of driving one node (labeled 0) are illustrated as a function of driving frequencies [(c)–(e)] and at four selected response nodes (labeled 1–4). Whereas at high driving frequencies responses are localized (algebraically in frequency, exponentially in internode distance), they are irregular across frequencies and among nodes in a resonance regime of intermediate driving frequencies. At low frequencies, globally homogeneous bulk oscillations emerge where the entire grid follows the driving signal (even if at only one node). All response strengths are quantified in terms of the amplitude of the frequency response relative to their low frequency limit as $ω \to 0$ $ω \to 0$ . From X. Zhang; cf. Zhang *et al.* (2019).

For large frequencies (compared to the range set by the Laplacian eigenvalues), the responses are strongly localized on the network, with amplitudes

c i (ω) \sim const \times ω - 2 d - 1

$c_{i} (ω) \sim const \times ω^{- 2 d - 1}$

(171)

in Eq. (159) asymptotically (as $ω \to \infty$ $ω \to \infty$ ) decaying exponentially with frequency and algebraically with graph theoretic distance $d = d (i, j)$ $d = d (i, j)$ response node $i$ $i$ has from the driven node $j$ $j$ . In Eq. (171) $const$ $const$ is a constant value, independent of frequency yet containing information about eigenvector overlaps and the identity of the driving and response nodes. This result complements the findings on localization reported in Fig. 31 for short-duration perturbations to signals with arbitrary frequency content and asymptotically quantifies localization. An exact asymptotic expansion up to order $ω - 2 d - 1$ $ω^{- 2 d - 1}$ shows that the entry $(i, j)$ $(i, j)$ of all powers $m < d (i, j)$ $m < d (i, j)$ of the Laplacian exactly equals 0, reflecting no existing paths from node $j$ $j$ to $i$ $i$ of length shorter than their topological distance $d (i, j)$ $d (i, j)$ ; see the Supplemental Material of Zhang et al. (2019).

For intermediate frequencies in and near the range defined by the Laplacian eigenvalues, resonances occur that induce spatiotemporal responses that are strongly inhomogeneous both as a function of frequency $ω$ $ω$ and among response nodes $i$ $i$ . Near resonance frequencies relative to response strengths substantially exceed those at low frequencies; see Fig. 34(b) as well as Fig. 2(f) of Zhang et al. (2019). Here the interaction network topology plays a major role in selecting which node and frequency responses are particularly strong (or weak), thereby stressing the distributed nature of these resonances. At low driving frequencies, the dynamics at all grid nodes follow the driving signal almost instantaneously, resulting in spatially homogeneous bulk oscillations.

Such an analysis transfers to fluctuating signals with distributed frequency content and distributed driven nodes, for instance, grid models driven by purely random processes as well as those driven with power frequency fluctuations characteristic of photovoltaic or wind power generators are well characterized by linear response theory. Only above 95% transmission line loads, a regime to be avoided in real operations for various reasons, does linear response theory yield substantial errors.

D. Blackouts as rare events due to large fluctuations

Dynamical instabilities of the power grid dynamics will become increasingly important as fluctuations from renewable resources become more frequent and possibly stronger. How close is a given system to an unstable state, and thus what is a typical timescale for fluctuations to kick the system out of its stable operating state? If strong fluctuations are relatively infrequent, it is desirable to quantify how rare they are. In addition, one would naively expect that fluctuations of power production with non-Gaussian distributions increase the risk of desynchronization, as large fluctuations are less suppressed in their broad tails. An analytic framework pursued by Hindes, Jacquod, and Schwartz (2019) quantifies the risk of escape depending on the deviations from Gaussian fluctuations. This work generalizes the analytic estimates of Schäfer et al. (2017) outlined in Sec. V.C.4 in several respects. It captures the distributed dynamics of large grids and Poissonian (rather than Gaussian) noise suited to model fluctuations in real-power production.

Hindes, Jacquod, and Schwartz (2019) employed a WKB approach applied to classical stochastic systems; cf. Dykman et al. (1994). They found that the rate of desynchronization may exponentially speed up or slow down, depending on how the statistics of fluctuations combines with the least stable mode of the network described by the Fiedler vector. In contrast to the Kramers-like formula for the escape rate [Eq. (107)], which depends only on the second moment, higher-order moments of the Poissonian noise are captured by the WKB framework, and their impact on the escape rate may be counterintuitive. As we later argue, the escape time $⟨ T ⟩$ $⟨ T ⟩$ is proportional to $exp [- (S (0) + Δ (n) S)]$ $\exp [- (S^{(0)} + Δ^{(n)} S)]$ , with the action $S$ $S$ expanded in powers of the distance to the bifurcation point. Higher-order corrections $Δ (n) S$ $Δ^{(n)} S$ to the leading term $S (0)$ $S^{(0)}$ can have both signs, increasing or decreasing the escape time. We summarize the assumptions for this WKB approach to sketch the main steps of deriving the escape time from a stable state. For detailed derivations see Hindes, Jacquod, and Schwartz (2019) and references therein.

Consider power grids described by the swing equation for synchronous machines, where an input power $¯ P i$ ${\bar{P}}_{i}$ is subject to fluctuations $p i (t)$ $p_{i} (t)$ . Given histogram data, such as from turbulence induced wind-generated power increments, we construct a time series $p i (t)$ $p_{i} (t)$ of fluctuations reproducing the increment histogram with bins of size $b$ $b$ . The ansatz is overdamped dynamics

˙ p i = - γ p i + ξ i (t),

${\dot{p}}_{i} = - γ p_{i} + ξ_{i} (t),$

(172)

with damping rate $γ$ $γ$ , and a statistical drive $ξ i (t)$ $ξ_{i} (t)$ , so large intermittent spikes are allowed, as they are observed for wind and solar sources. The stochastic drive is given as $ξ i (t) = \sum b n g i b δ (t - t i b [n])$ $ξ_{i} (t) = \sum_{b n} g_{i b} δ (t - t_{i b} [n])$ , where the power increment $g i b$ $g_{i b}$ at $i$ $i$ is the average of the pulse amplitude over the bin $b \in {1, 2, \dots, M}$ $b \in {1, 2, \dots, M}$ of the histogram, where averaged amplitudes are the measured increments $g i = p i (t + 1) - p i (t)$ $g_{i} = p_{i} (t + 1) - p_{i} (t)$ and $t i b [n]$ $t_{i b} [n]$ denotes the time at which the $n$ $n$ th such increment in bin $b$ $b$ at unit $i$ $i$ occurs. For modeling, the driving signal is assumed to be Poisson shot noise, so the time between two events where unit $i$ $i$ receives a power increment within bin $b$ $b$ is exponentially distributed with a rate $ν i b$ $ν_{i b}$ that is consistent with the histogram.

Consider a synchronized phase-locked state that emerges through an inverse saddle-node bifurcation as the coupling strengths increase such that unstable, saddle phase-locked states are in the vicinity of it. Moreover, take typical fluctuations to be small compared to distances to the saddles and large fluctuations to be rare. Such fluctuations are not captured in the large-deviation approach of Nesti, Zocca, and Zwart (2018) described in Sec. VII.B.3, where an overload may also be due to the many small deviations that accumulate.

The thinking behind the analytical approach of Hindes, Jacquod, and Schwartz (2019) is that a large fluctuation drives the system along a most-likely path in phase space. The path connects the stable phase-locked state with the saddle, while fluctuations and the network dynamics coact as to maximize the probability of desynchronization. To derive the optimal path from classical mechanics, a generalized Fokker-Planck equation for the stochastic network dynamics is adapted to incorporate Poissonian noise. It determines the probability distribution $ρ$ $ρ$ of finding oscillator phases $Φ$ $Φ$ , phase velocities $v$ $v$ , and fluctuations $p$ $p$ at time $t$ $t$ . To analyze large fluctuations that are rare, we project on the exponential tails of $ρ$ $ρ$ and insert a WKB ansatz

ρ (Φ, v, p, t) = B exp {- S (Φ, v, p, t)}

$ρ (Φ, v, p, t) = B \exp {- S (Φ, v, p, t)}$

(173)

into the Fokker-Planck equation. Next the action $S (X + δ X)$ $S (X + δ X)$ with $X \equiv (Φ, v, p)$ $X \equiv (Φ, v, p)$ is expanded in deviations $δ X$ $δ X$ about $X$ $X$ (here the stationary state), where $S (X) ≫ 1$ $S (X) ≫ 1$ , while $| δ X | \propto 1$ $| δ X | \propto 1$ is assumed and only leading-order terms in the first derivatives $\partial Φ S$ $\partial_{Φ} S$ , $\partial v S$ $\partial_{v} S$ , and $\partial p S$ $\partial_{p} S$ are retained. As a remark, this ansatz allows one to analyze the shape of the probability distribution even in the remote tails for which $δ X \propto Ω ≫ Ω 1 / 2$ $δ X \propto Ω ≫ Ω^{1 / 2}$ if $Ω 1 / 2$ $Ω^{1 / 2}$ is the typical size of a fluctuation in the stationary state, in contrast to the typical and small fluctuations considered in Eq. (106) in Kramers’s escape theory or entering Eq. (150) as it is derived from large-deviation theory. The insertion of Eq. (173) into the Fokker-Planck equation yields a Hamilton-Jacobi equation for an “action” $S (Φ, v, p, t)$ $S (Φ, v, p, t)$ , from which a classical Hamiltonian $H (Φ, v, p, λ Φ, λ v, λ p)$ $H (Φ, v, p, λ^{Φ}, λ^{v}, λ^{p})$ can be read off. The Hamiltonian depends then on “coordinates” $Φ$ $Φ$ , $v$ $v$ , and $p$ $p$ , and their conjugate momenta $λ Φ, λ v, λ p$ $λ^{Φ}, λ^{v}, λ^{p}$ equates to zero [ $H (Φ, v, p, λ Φ, λ v, λ p) = 0$ $H (Φ, v, p, λ^{Φ}, λ^{v}, λ^{p}) = 0$ ]. The action reduces to

S (Φ, v, p) = \sum i [\int λ Φ i d Φ i + \int λ v i d v i + \int λ p i d p i] .

$S (Φ, v, p) = \sum_{i} [\int λ_{i}^{Φ} d Φ_{i} + \int λ_{i}^{v} d v_{i} + \int λ_{i}^{p} d p_{i}] .$

(174)

The next goal is to determine the optimal (least action) path in phase space as a solution to the classical Hamilton’s equations of motion, subject to the boundary conditions to connect the stable fixed point with the unstable saddle. The action along this optimal path is stationary [therefore, we skipped the explicit time dependence of $S$ $S$ in Eq. (174)] and enables us to estimate the expected waiting time $⟨ T ⟩$ $⟨ T ⟩$ for desynchronization via

ln ⟨ T ⟩ \sim S (Φ s, 0, 0) + const,

$\ln ⟨ T ⟩ \sim S (Φ^{s}, 0, 0) + const,$

(175)

where $Φ s$ $Φ^{s}$ denotes the stationary solution. To estimate $S (Φ s, 0, 0)$ $S (Φ^{s}, 0, 0)$ , the solution is determined for coupling strengths $K = K SN (1 + κ)$ $K = K_{SN} (1 + κ)$ close to $K SN$ $K_{SN}$ , where saddle-node bifurcation happens. To lowest order in $κ$ $κ$ , parametrizing the distance to the bifurcation point, the solution on the one-dimensional submanifold takes the form

Φ i (t) = Φ SN i + C κ 1 / 2 r i x (t),

$Φ_{i} (t) = Φ_{i}^{SN} + C κ^{1 / 2} r_{i} x (t),$

(176)

with $x \in [- 1, + 1]$ $x \in [- 1, + 1]$ such that, for $x = - 1$ $x = - 1$ , $Φ i (t) = Φ SN i - C κ 1 / 2 r i = Φ * i$ $Φ_{i} (t) = Φ_{i}^{SN} - C κ^{1 / 2} r_{i} = Φ_{i}^{*}$ . The stable fixed point is the starting point, and for $x = + 1$ $x = + 1$ we have $Φ i (t) = Φ SN i + C κ 1 / 2 r i = Φ s i$ $Φ_{i} (t) = Φ_{i}^{SN} + C κ^{1 / 2} r_{i} = Φ_{i}^{s}$ as the unstable saddle as the end point of the optimal path; $r i$ $r_{i}$ is the component $i$ $i$ of the Fiedler mode of the Laplacian that encodes the network topology. This mode is most sensitive to external fluctuations and mediates the unlikely escape. The constant $C$ $C$ is independent of $x i$ $x_{i}$ and $κ$ $κ$ and depends on the grid adjacency matrix, the phases at the bifurcation point $Φ SN i$ $Φ_{i}^{SN}$ , and the Fiedler mode $r i$ $r_{i}$ .

As it turns out (Hindes, Jacquod, and Schwartz, 2019), to lowest order in $κ$ $κ$ the action $S (Φ s, 0, 0)$ $S (Φ^{s}, 0, 0)$ depends on the damping constant $γ$ $γ$ , the saddle-node coupling $K SN$ $K_{SN}$ , and the damping rate of power fluctuations squared $(γ 2)$ $(γ^{2})$ . Moreover, each node of the grid $i$ $i$ contributes to $S$ $S$ , a term proportional to $μ - 1 2$ $μ_{2}^{- 1}$ , with $μ 2 = \sum b ν i b g 2 i b$ $μ_{2} = \sum_{b} ν_{i b} g_{i b}^{2}$ the fluctuation variance. This means that to lowest order in the distance $κ$ $κ$ from the saddle, the expected time to blackout $⟨ T ⟩ \propto e - S$ $⟨ T ⟩ \propto e^{- S}$ is insensitive to higher-order moments $μ n$ $μ_{n}$ , $n \geq 3$ $n \geq 3$ , of the fluctuations (for Poissonian noise). However, at higher orders in $κ$ $κ$ , the corrections to $S$ $S$ and $Δ (n) S$ $Δ^{(n)} S$ scale with $κ$ $κ$ according to $κ n - 1 / 2$ $κ^{n - 1 / 2}$ . They depend on higher-order moments $μ n$ $μ_{n}$ and on the balance between the positive and negative Fiedler’s mode components. For $n = 3$ $n = 3$ , it is the product of the skewness of the power fluctuation distribution with the skewness of the Fiedler mode that determines whether the non-Gaussianity leads to an increase or decrease of the desynchronization rates captured by $Δ (n) S$ $Δ^{(n)} S$ . Symmetric power fluctuations with $μ 3$ $μ_{3}$ close to zero yield an increase in the desynchronization rate, as naively expected. Overall, this result (Hindes, Jacquod, and Schwartz, 2019) demonstrates how the network topology interacts with the stochastic dynamics in a nontrivial way.

A further application of the relation between desynchronization and the Fiedler-mode values amounts to a dimensional reduction of the phase space. The phase space is usually high dimensional for a power grid, but the possible desynchronization paths lie in a low-dimensional subspace under the given assumptions. For networks in which the saddle-node bifurcation is induced by a single overloaded edge with phase difference $π / 2$ $π / 2$ [see also Fliscounakis et al. (2013), Manik et al. (2014), and Rohden et al. (2017)], a so-called synchronized subgraph approximation becomes exact (Hindes, Jacquod, and Schwartz, 2019). The Hamilton’s equations of the two subsystems that desynchronize reduce to a single noisy oscillator system in relative phase-space coordinates. More generally, the subgraph approximation is sensible if the corresponding network partition is guided by approximately uniform Fiedler-mode values. Thus, rare desynchronization events become analytically predictable in spite of the high-dimensional phase space of the grid dynamics.

E. Ensembles of power grids for distributed fluctuations

In Secs. VIII.A–VIII.C the focus was on the spatiotemporal propagation of a perturbation applied to a single or a few nodes, or the addition or removal of individual lines. Mureddu et al. (2015) introduced a suitable framework for estimating the overall energy mismatch between day-ahead estimates and real data of an entire grid that should be balanced with an appropriate trade on the energy balancing market. Differences between production and consumption result from perturbations at all nodes simultaneously. The framework was later employed in different applications (Mureddu and Damiano, 2017; Korjani et al., 2018; Mureddu and Meyer-Ortmanns, 2018).

1. Basic approach

The idea is to consider an ensemble of power grids in analogy to an ensemble of microstates in statistical mechanics. The power grids differ in their “microstates,” which enter average values of global observables. Global observables may be the market volume, the amount of energy that the energy balancing market must compensate for, the costs for this amount of energy (Mureddu and Meyer-Ortmanns, 2018), the resilience of the grid (Mureddu and Damiano, 2017), or the optimal positioning of storage devices (Korjani et al., 2018).

The individual “configurations” that represent the microstates are generated from a first reference configuration of producers and consumers that is representative of part of the considered region, a certain partition between renewable and conventional power generation, a certain time during the day, or a season of the year. The reference configuration is chosen to satisfy the OPF equations; it is regarded as the operation point of the system. For either transmission or distribution grids, the ensemble of microstates is generated by applying fluctuations in production and consumption relative to the reference configuration.

The fluctuations can be due to load fluctuations, forecast errors for renewables, intraday electricity trading, or other factors. In the simplest case, the fluctuations are chosen from a truncated Gaussian distribution. Gaussian-like forecast errors are defined by standard deviations $σ i$ $σ_{i}$ at nodes $i$ $i$ and represent the expected power variations at that node at a given time. For a load of type $ℓ$ $ℓ$ (where $ℓ \in {w, pv, \dots}$ $ℓ \in {w, pv, \dots}$ denotes wind, photovoltaic, and other types of power generation) and power demand $P ℓ$ $P_{ℓ}$ , the standard deviation of the forecasting error is denoted as $σ ℓ$ $σ_{ℓ}$ , and $m ℓ$ $m_{ℓ}$ and $M ℓ$ $M_{ℓ}$ are the minimum and maximum values of the support of the distribution $ρ ℓ$ $ρ_{ℓ}$ of $P ℓ$ $P_{ℓ}$ , respectively, representing power constraints of the different generators. Depending on the different assignments $P ℓ$ $P_{ℓ}$ to nodes $i$ $i$ , a microstate is sampled by adding a random value to the expected power production or consumption at every node $i$ $i$ , where the random variable is extracted from the truncated probability density function $ρ t i (x)$ $ρ_{i}^{t} (x)$ as follows:

ρ t i (x) = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 0 if x < m i, ρ i (x) if m i < x < M i, 0 if x > M i,

$ρ_{i}^{t} (x) = {\begin{cases} 0 & if x < m_{i}, \\ ρ_{i} (x) & if m_{i} < x < M_{i}, \\ 0 & if x > M_{i}, \end{cases}$

(177)

where the original densities (with support in all of $R)$ $R)$ are $ρ i (x) = (2 π σ 2 i) - 1 / 2 exp (- x 2 / 2 σ 2 i)$ $ρ_{i} (x) = {(2 π σ_{i}^{2})}^{- 1 / 2} \exp (- x^{2} / 2 σ_{i}^{2})$ for wind, photovoltaic, or other generation at nodes $i$ $i$ . For skewed distributions, modeling, e.g., power production by wind, we may choose the following Weibull distribution:

ρ Weibull (x; λ, a) = {a λ (x λ) a - 1 e - (x / λ) a if x \geq 0, 0 otherwise,

$ρ^{Weibull} (x; λ, a) = {\begin{cases} \frac{a}{λ} (\frac{x}{λ})^{a - 1} e^{- (x / λ)^{a}} & if x \geq 0, \\ 0 & otherwise, \end{cases}$

(178)

with $a = 2$ $a = 2$ (Lun and Lam, 2000; Seguro and Lambert, 2000; Lu et al., 2013) and $λ = 2 P w / \sqrt π$ $λ = 2 P_{w} / \sqrt{π}$ , where $P w$ $P_{w}$ represents the wind power production of the reference configuration and $Γ$ $Γ$ denotes the gamma function. Specific values of the parameters fixing the respective distribution depend on the load and the type of renewable energy and are chosen from recorded data.

As a next step, quantities like the resulting mismatch $P j i$ $P_{i}^{j}$ of the power at node $i$ $i$ in configuration $j$ $j$ are measured and the $P j i$ $P_{i}^{j}$ summed over all nodes. This way we obtain the total mismatch $P j$ $P^{j}$ in power production (the so-called market volume) for a given configuration $j$ $j$ , which the market will balance. To obtain representative values of the market volume, we sample a sufficiently large number of configurations to include in the ensemble. In particular, one can analyze how the distribution of the market volume over the volume size depends on the distribution of fluctuations (normal or Weibull), and (via the choice of the reference configuration) on the time of day, the season, the geographical zone, and the percentage of renewables with respect to the total production (Mureddu and Meyer-Ortmanns, 2018). The distribution of this market volume enters the market costs, the average prices, and the average profit per technology that is involved in the production.

2. A second network layer

To model the energy market and to determine market costs, prices, and profit, a second network layer describes market trading. Here one option is agent-based modeling with network nodes representing agents [cf. Han et al. (2019)] such that the agents are the retailers, with one retailer per conventional power station of the physical grid. In the work of Mureddu and Meyer-Ortmanns (2018), the agents first have to undergo a learning phase, in which they learn how to place optimal bids in competitive auctions with the aim of buying (or selling) in the most profitable way. To simulate how real market operators acquire knowledge about the market in time and adapt their decisions, a modified Roth-Erev algorithm (Nicolaisen, Petrov, and Tesfatsion, 2001; Mureddu et al., 2015) adjusts the offer propensities of agents in a self-consistent way, with the goal of maximizing profits. The agents interact via a so-called market authority that provides a link between the physical grid and the market. The market authority knows the mismatch from the physical grid and takes the bids of the retailers, accepts or rejects these bids, informs the retailers about the decision, and proceeds until the required mismatch in energy is covered at the lowest possible costs.

If the energy balance is restored with energy provided by the subset of retailers who offered the energy at the lowest price, the feedback on the physical grid stability must be checked since the economically best selection need not be reasonable from the viewpoint of grid stability. Lines might get overloaded if the conventional generators that were selected for selling the energy happen to spatially cluster together. Thus, this framework of different network layers, coupled via a market authority, allows one to analyze the feedback from economical (low costs) to physical (high stability) optimization objectives.

The physical grid stability is particularly endangered if the retailers behave like arbitrageurs when the reserve energy price falls below the price on the intraday market. In such a situation, retailers play a kind of minority game (Ritmeester and Meyer-Ortmanns, 2021). Modeling their behavior accordingly leads to suggestions on how to control arbitrage. For example, if the few large parties contributing to the market (rather than the many small ones) are made risk averse due to small penalties, it has a disproportionately large effect on reducing the abuse of price differences in terms of arbitrage. Moreover, from the remarkable analogy of the minority game with spin glasses it becomes understandable why a larger number of retailers may not at all reduce the fluctuations in arbitrage, as one would naively expect. Related features of underlying phase transitions are also visible in realistic markets.

3. Alternatives for treating uncertainties

A number of alternative approaches exist to deal with inherent uncertainties in power and demand. The goal is to derive distributions of induced fluctuations. Induced fluctuations may refer to the energy mismatch (as previously discussed), or to induced voltage or frequency fluctuations. An approach called chance-constrained ac optimal power flow (ac CC OPF) was described by Roald and Andersson (2018) and references therein. In contrast to the optimal power flow summarized in Sec. III.B, chance constraints ensure that the system constraints will be satisfied only with a specified probability. The framework of ac CC OPF optimizes not only the scheduled dispatch with respect to costs but also the procurement of reserve and voltage control during deviations from expected values. As outlined by Roald and Andersson (2018), a number of different options exist for solving this optimization problem under uncertainty. The first one is a one-shot optimization, where the optimal solution is found while respecting all constraints simultaneously (which is demanding in view of the complexity of the problem). The second one is an iterative solution algorithm. The third option similar to the previously mentioned ensemble approach employs Monte Carlo simulations for deriving uncertainty margins and samples uncertainties in power and demand. Here the resulting power flows are calculated for a large number of sample realizations drawn from a given ensemble, leading to Monte Carlo–based uncertainty margins for the voltage. For a comparison of the efficiency of the various approaches, see Roald and Andersson (2018).

IX. Summary and Outlook

Modern power grids are tasked with incorporating increasing shares of renewable energy sources. This ongoing drastic transformation is crucial to mitigating the climate crisis but poses novel challenges to our understanding of the collective dynamics of power grids. Although the underlying equations have been known for a long time, the treatment of such grids on the systems level remains difficult, mainly due to their nonlinearities, distinct heterogeneities in the dynamics of the nodes, and their interaction topology, as well as the various kinds of perturbations and fluctuations present. Additionally, power grids are embedded in and interconnected to several other complex systems, such as energy markets, consumers, and weather. We discussed here central challenges of decentralized and heterogeneous power grids and ways to analyze, understand, design, and influence them.

The review started with a description of power grids based on the principles of physics and the treatment of the most important node models. A crucial aspect here is that renewable sources are connected via electronic inverters. This necessitates an extension of the classic models and leads to fundamental changes in the dynamics of power grids. To study these models and the variety of analysis techniques, datasets on dynamical parameters and network topologies are necessary. However, data on real-world power grids are scarce and often incomplete. Therefore, we presented in Sec. IV typical classes of synthetic models that mimic the main properties of real grids, such as the classic IEEE test cases, but also more recent developments on generating novel network ensembles, which are often closer to real grids and emphasize the physics perspective. Section V discussed an elementary grid containing only one transmission line. This simple bistable system served to introduce the most important static solutions and voltage limits, the bifurcations of the dynamics, and recent stability concepts such as basin stability and stochastic stability, providing a deeper understanding of stability with respect to large and sustained perturbations.

In Sec. VI, recent achievements and insights for realistic networks were presented. To function as a power grid, the networked dynamics needs to have a stable synchronous (phase-locked) state. However, power grids are typically multistable, exhibiting multiple synchronous yet also asynchronous states that are unsuitable for grid operation. This multistability increases substantially if the ubiquitous losses are included, leading to new types of both synchronized and desynchronized states. We gave an overview of the most successful analytic and probabilistic approaches for managing such a rich variety of dynamics, as well as classic linearization-based methods.

The structural stability of grids was the topic of Sec. VII. In complex grids the outage of a single transmission element may induce cascades leading to a large-scale blackout. Such events are well documented in real grids all over the globe, among them several “monster blackouts.” When a single element fails or is added, the network flows reroute in response. A main grid parameter influencing this rerouting is the network distance. But the splitting of a grid into different communities also has a strong impact. New mathematical descriptions of the rerouting process combining methods from statistical physics, nonlinear dynamics, and graph theory are capable of describing these reroutings in new ways and reliably predicting critical links in a network. It turns out that a local stability analysis is often not sufficient and may even be misleading. The occurrence of cascades following individual failures can also be enhanced by transient effects, which shows that dynamical and structural aspects are deeply interwoven. Furthermore, reroutings may yield counterintuitive consequences if transmission elements are added or removed. For instance, under certain conditions the reinforcement of a transition line or the addition of a new line may induce a loss of capacity and stability, a phenomenon known as Braess’s paradox.

The power grid is subject to fluctuations. Section VIII discussed this essential problem with a particular focus on the characteristics introduced by renewables. The non-Gaussian nature of their fluctuations, often characterized by heavy-tailed distributions, such as those due to the cloud structure affecting photovoltaic generation, calls for new modeling and analysis approaches. First, a linear response theory was developed to describe the collective spatiotemporal grid dynamics subject to both stationary and transient and nonstationary and distributed input fluctuations. This led to an efficient way of identifying vulnerability patterns. Second, the emergence of resonances, bulk oscillations, or localizations was explained. Recent research demonstrated using a WKB ansatz that even desynchronization events in grids, which originate from strong but rare fluctuations, can be predicted. The results revise the naive expectation that fat-tail distributions in fluctuating power production always increase the number of rare events of desynchronization. In real networks with renewable components, fluctuations in generation and consumption occur in parallel; i.e., one cannot restrict the robustness analysis to the failure of a single line. Statistical physics provides alternative approaches to treating this problem, such as the consideration of an entire ensemble of power grids, differing in the realization of fluctuations in its microstates, or by chance-constrained ac optimal power flow.

Aspects of power system control and monitoring were only mentioned in this review but will become increasingly important in future power grids with many distributed, fluctuating power sources and low inertia (Milano et al., 2018). This field of research spans from single machines to the entire energy system. Current challenges are the development of control systems for grid-forming inverters that guarantee global dynamic stability (Colombino et al., 2019; Schiffer, Efimov, and Ortega, 2019) or virtual inertia systems that may replace the mechanic inertia of synchronous machines (Chen et al., 2011; Kerdphol et al., 2019). On the grid level, researchers thrive for a better understanding of the interplay of control systems and network dynamics (Tumash, Olmi, and Schöll, 2019; Totz, Olmi, and Schöll, 2020) and the development of new concepts for the control of complex networks (Cornelius, Kath, and Motter, 2013; Huang et al., 2019). Finally, the operation of the large-scale load-frequency control system is subject to multiple external influences, including markets and regulations (Kruse, Schäfer, and Witthaut, 2021a, 2021b).

This review demonstrated that the analysis and design of decentralized power grids demand an interdisciplinary approach where techniques and concepts of power engineering and control theory are newly combined with those from statistical physics, complex dynamical systems, and network science. Such an approach will continue to form the backbone for gaining a deeper understanding of all aspects discussed here and will provide the inspiration for important future research directions.

Future integrated studies will have to connect the power grid to a broader range of important aspects of the energy system, such as the energy market, the information and communication infrastructure driving the control, consumer behavior, electromobility, political and planning constraints, etc. All these aspects point toward the nature of the power grid as part of a multilayered network that is embedded in and connected to other complex systems. It is becoming increasingly clear that machine learning techniques have enormous potential in the study and control of complex heterogeneous systems. Specifically, graph neural networks have vast untapped potential. The question of how to effectively combine such data-driven approaches with physics, engineering-based concepts, and generic modeling remains open in the context of power grids.

Finally, a deeper understanding of fundamental concepts underlying the collective dynamics of power grids should enable more than an energy transition in highly developed countries with well-connected grids required to stop the worsening of the climate crisis. It may also enable an implementation of sustainable power systems in rural areas of the global South and in megacities that now may be designed from scratch. Finally, it should improve an understanding of small-scale energy islands that can play a crucial role in the sustainable electrification of local communities that still remain off grid today.

List of Symbols and Abbreviations

CC OPF

chance-constrained optimal power flow

DC OPF

linearized or dc optimal power flow

EMF

electromotive force

ENTSO-E

European Network of Transmission System Operators for Electricity

HVDC

high voltage directed current

LODF

line outage distribution factor

OPF

optimal power flow

PI law

proportional-integral control law

PTDF

power transfer distribution factor

pu (system)

per unit (system of rescaled units)

RoCoF

rate of change of frequency

WKB

Wentzel-Kramers-Brillouin

Acknowledgments

This review was produced during our work on the collaborative research project “Kollektive Nichtlineare Dynamik Komplexer Stromnetze” (CoNDyNet), which was funded by the German Federal Ministry of Education and Research (BMBF Grants No. 03EK3055 and No. 03SF0472). D. W. acknowledges support from the Helmholtz Association via the joint initiative “Energy System 2050—A Contribution of the Research Field Energy” and Grant No. VH-NG-1025. J. K. was supported by Russian Ministry of Science and Education Agreement No. 075-15-2020-808.

Reviews of Modern Physics

Collective nonlinear dynamics and self-organization in decentralized power grids

Dirk Witthaut, Frank Hellmann, Jürgen Kurths, Stefan Kettemann, Hildegard Meyer-Ortmanns, and Marc Timme

Rev. Mod. Phys. 94, 015005 – Published 28 February 2022

Abstract

Physics Subject Headings (PhySH)

Authors & Affiliations

Article Text